— 


a c Further Considerations on the Statistical Mechanics of Biological 


THE BULLETIN OF 
Be ctical’. 


BIOPHYSICS 


JUNE 1959 


A Relational Theory of Biological Systems II—Robert Rosen - - - 109 


A Contribution to the Study of Diffusion of Neutral Particles Through 


Penne Petras Pen Cie. hse ee a) ap afi a a eet a 129 
A Note on the Transient Stage of the Random Dispersal of Logistic 

Populations——Richard Barakat - ----++-+-+-++-+-- 141 
A Note on Population Growth Under Random Dispersal—H. D. 

ST ee ee 153 
Some Remarks on the Mathematical Theory of Nutrition of Fishes— 

N. Rashevsky- - ---+----++ +++ eee eee eee 161 
A Note on the Nature and Origin of Life——-N. Rashevsky - - - - - 185 


Feedback in Physiological Systems: An Application of Feedback 
Analysis and Stochastic Models to Neurophysiology—Alan R. 


Wdolpheareen nie msm win ba oe 2 6 2 shee ee ee 195 


Associations Edward H. Kerner - - - +--+ ++ eee 22% 


THE UNIVERSITY OF CHICAGO PRESS - CHICAGO 


VOLUME 21 Brag res. auth NUMBER 2 


Mf 


is 


The Bulletin of 
MATHEMATICAL BIOPHYSICS 


Ep1tor: 
N. RASHEVSKY, UNIVERSITY OF CHICAGO 
AssoctaTE Epirors: 
H. D. LANDAHL, UNIVERSITY OF CHICAGO 
ANATOL RAPOPORT, UNIVERSITY OF MICHIGAN 


The BULLETIN is devoted to publications of research in Mathematical 
Biology, as described on the inside back cover. 


AUTHORS: 


Please read carefully the inside back cover. 


Published March, June, September, and December. Subscription 
rates per volume: U.S.A., $10.00; Canada and Pan America, $10.50; 
all other countries in Postal Union, $11.00. Single copies: $3.25. 
Checks should be made payable to The University of Chicago Press, 
5750 Ellis Avenue, Chicago 37, Illinois. Subscriptions will be en- 
tered to start with the first issue to be published after receipt of 
order. Editorial correspondence: The Editor, The Bulletin of Math- 
ematical Biophysics, 5741 Drexel Avenue, Chicago 37, Illinois. 
Subscribers are requested to notify the Press and their local post- 


master immediately of change of address, giving both the old and 
the new address. 


Second-class postage paid at Lancaster, Pennsylvania. 


[Copyright 1959 by the University of Chicago] 


BULLETIN OF 
MATHEMATICAL BIOPHYSICS 
VOLUME 21, 1959 


A RELATIONAL THEORY OF BIOLOGICAL SYSTEMS II 


ROBERT ROSEN 
COMMITTEE ON MATHEMATICAL BIOLOGY 
THE UNIVERSITY OF CHICAGO 


The general Theory of Categories is applied to the study of the (M, R)- 
systems previously defined. A set of axioms is provided which character- 
ize ‘‘abstract (M, R)-systems,’’ defined in terms of the Theory of Categories. 
It is shown that the replication of the repair components of these systems 
may be accounted for in a natural way within this framework, thereby 
obviating the need for an ad hoc postulation of a replication mechanism. 

A time-lag structure is introduced into these abstract (M, R)-systems. 
In order to apply this structure to a discussion of the ‘‘morphology’’ of 
these systems, it is necessary to make certain assumptions which relate 
the morphology to the time lags. By so doing, a system of abstract biol- 
ogy is in effect constructed. In particular, a formulation of a general 
Principle of Optimal Design is proposed for these systems. It is shown 
under what conditions the repair mechanism of the system will be local- 
ized into a spherical region, suggestive of the nuclear arrangements in 
cells. The possibility of placing an abstract (M, R)-system into optimal 
form in more than one way is then investigated, and a necessary and 
sufficient condition for this occurrence is obtained. Some further impli- 
cations of the above assumptions are then discussed. 


In a previous paper (Rosen, 1958b, hereinafter referred to as (II)) 
we outlined a possible approach to the theory of general systems, 
based on the notion of a category. The present work is an attempt 
to indicate how such an approach may, through the theory of (M, R)- 
systems introduced by us (Rosen, 1958a, hereafter referred to as 
(I)), be made to bear more directly on a number of important bio- 
logical problems. We shall assume that the reader is familiar with 
the material contained in (I) and (II), and we shall use the notation 
and terminology of those papers without further comment. 


I. Abstract (M, R)-systems. Our first task is to formalize the no- 
tion of (M, R)-system as described in (I), in terms of the more ab- 
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stract theory which we developed in (II). We recall that an (M, R)- 
system consists essentially of a system M, together with certain 
auxiliary systems, designated by &,, such that each component 
M, of M has associated with it exactly one of the #;, such that the 
output of R, is a replica of the component M, with which it is as- 
sociated. The inputs to each F#, consist of a subset of the set of 
environmental outputs of M. We also imposed certain other struc- 
ture on the (M, R)-systems in (I), in connection with the lags of the 
system; we shall deal with these properties in the next section. 
For the moment, however, we shall ignore the time lag structure, 
and consider an (M,R)-system as an assembly of components M, 
with their associated systems Ff, as described above. 

It may be useful to have a concrete example before our eyes, on 
which we may explicitly illustrate the techniques which we shall 
develop below. To this end, let us consider the specific (M, R)- 
system which is labeled Figure 1 below (this is the same diagram 
as the one displayed on p. 254 of (I)); this figure represents the 
ordinary block diagram of the given system. By our discussion in 
the Introduction of (II), we know that this diagram is inherently 
ambiguous, in the sense that, when two arrows are shown as origi- 
nating from a component (such as M, of the figure) we cannot tell 
whether these arrows represent different outputs of the component, 
or merely the same output being supplied to two different compo- 
nents. We shall henceforth assume that each component in this 
diagram which emits two arrows actually produces two distinct 
outputs. 


FIGURE 1 
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Let us now return to the general discussion. By Theorem 2 of 
(II) we know that we can find an abstract block diagram consisting 
of a collection of sets and mappings from a suitable subcategory 
of the category S of sets, which will represent a given system M. 
Further, we may suppose that this abstract block diagram is in the 
canonical form described by Theorem 4 of (II); in this case, each 
component of the system M under consideration is represented by a 
single mapping, the domain of which is contained in a cartesian 
product of sets of S. It will be found that this particular canonical 
representation will be most convenient both for the construction of 
our representation of abstract (M,R)-systems and for our subse- 
quent discussion of the time-lag structure of these systems. 

To see explicitly how to construct such canonical abstract block 
diagrams, we carry out this construction in detail for the special 
example mentioned above. First, we observe that, once the am- 
biguity concerning multiple outputs is resolved, we may apply 
Theorem 4 of (II) directly to ordinary block diagrams; we thereby 
obtain from a given ordinary block diagram a new ordinary block 
diagram in which each component emits a single output. Figure 2 
below shows the diagram which results when this process is ap- 
plied to the set of components M, of Figure 1; the notation is 
chosen in the obvious manner. The emission of two arrows from a 
single component in this canonical diagram will now be seen to 
always have an unambiguous meaning; namely, the same output is 
to be provided to two components. 


FIGURE 2 
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In Figure 3 below, we show the abstract block diagram obtained 
from the diagram of Figure 2 above by the procedure developed in 
(II). The notation is chosen so that the mapping corresponds to 
the component M,* of the ordinary diagram; the set A ,,n represents 
the output of the component M_ which serves as an input to the 
component M, (regardless of the superscripts of these components). 
Sets with only a single subscript represent environmental inputs 
and outputs. It will be seen that this abstract diagram is already 
in the canonical form of Theorem 4 of (II). 


A A 
e 5 
A, A3 A\ jk A, A; 
NA ee 
23 


A % t4 "i 
ae f t! A 
1 —~ z8 
ee ogee =, eps 
1 
\ fk \ jk 
As Axo 
FIGURE 38 


We must now find a means of representing the systems R, in 
terms of the Theory of Categories. From our previous discussions 
in (J) and (II), it follows that the &, must be represented by map- 
pings of some type. An inspection of the definition of (M, R)- 
system shows that the domain of these mappings must be repre- 
sented by a cartesian product of sets which are elements of the 
set © of environmental outputs of the system M (it will be shown 
below that @ is a set of sets), and the ranges of these mappings 
must contain the mappings which represent the related component M,. 

In greater detail, let us suppose that M is represented by an ab- 
stract block diagram in the manner previously described. We ob- 
serve at the outset that certain sets in the abstract block diagram 
of M are expressible as cartesian products of other sets in the 
diagram; such sets will be referred to as decomposable. Sets of 
the diagram which cannot be expressed as a cartesian product of 
other sets in the diagram will be called indecomposable. For ex- 
ample, the set which is the domain of the mapping f, in Figure 3 
above is decomposable; all other sets in this diagram are inde- 
composable. In particular, the canonical form which we have im- 
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posed on the diagram of M implies that the range of every mapping 
of the diagram is indecomposable (although the domains of these 
mappings are in general decomposable). In fact, it is not difficult 
to prove that this last is a necessary and sufficient condition for 
an abstract block diagram to be in the canonical form under 
discussion. 

Let the indecomposable sets in the abstract block diagram of M 
be enumerated in some fashion; say as A, aM Ser A,. Then 
every other set in the diagram will, according to Postulate A.B.D. 


2 of (Il), necessarily be of the form J] Ay where & runs through 
keK 


some finite index set K, and A, is indecomposable for each keK. 
If now f is a mapping in the abstract block diagram of M, the do- 
main and range of which we shall denote respectively as d(f) and 
r(f), then by Postulate A.B.D. 1 of (Il) there exists a set A= 


IT A, such that d(f) C A, and also an indecomposable set A; 
keK 
such that r(f) ¢ A; We shall find it convenient in the sequel to 


regard f as a mapping in the set of mappings H(A, A;.) (cf. (1), 
p. 321), even though f is not defined on all of A in general; we may 
accomplish this in a rigorous manner if we adopt the following 
artifice (which will be more thoroughly discussed in the next sec- 
tion): let f’ be a mapping in H(A, A; ) with the property that on the 
set d(f), f =f’. We shall require that for each such mapping f’ the 
following condition shall be fulfilled: if 2 is any element of the set 
A -d(f), then f’(z) has an infinite operation lag. The usefulness 
of this artifice will become apparent as we proceed with our repre- 
sentation. 

In (1), we denoted the collection of environmental outputs of the 
system M as ®. According to the representation theory developed 
in (II), @ will be represented by a collection of (indecomposable) 
sets in the abstract block diagram of M. Thus in Figure 3 above, 
® consists of the sets A,, A,, 43, 44; As, A,, A,, A,- Let the 
sets which comprise @ be enumerated in some fashion as BAG 
Bey iee; Db, Ther tie totality of possible inputs to the mappings 
which are to represent the components F, of an (M, R)-system con- 
sist of the elements of the sets of the form IT B; where J is a 


jeJ 
suitable index set, and B, «@ for each jeJ. To make this clear, let 
i ae: 
us consider once again the (M, R)-system shown in Figure 1. The 
input to the component F,,7=1, 2,..., 8; of this system is there 
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denoted by @,; according to our representation theory, we find for 
example that the set ©, = {p,, p,, pg} is to be represented in the 
corresponding abstract block diagram by the cartesian product 
A, x A, x A, of Figure 3 above; likewise the other sets @, are 
similarly represented by cartesian products in the more abstract 
situation. 

We may now undertake the definition of the representatives of 
the R; components. If M, is a component of a system M such that, 
in the abstract block diagram of M the component M, corresponds to 
a mapping f, where f:A — A; then we assign to the associated 
R ,-system a mapping 


®,: I] B;,— H(A, 4;,). 
jeJ 


In Figure 3 above, for example, this means that corresponding to 
the mapping f,.9 say, we have a mapping 


© ovASer A; x AL H(A, ,,4,) 


fg,o 


and corresponding to the mapping f, we have a mapping 


6,4? 4,5)° 


O, :4, H(A, xA A 


The reader may find it helpful to construct explicitly the mappings 
®, corresponding to the remaining mappings f of Figure 3 above, 
and even to construct similar examples of his own, in order to 
better familiarize himself with the mechanics of this construction. 

We thus find that we may construct an abstract representation of 
an (M, R)-system in the following manner: we construct the abstract 
block diagram (in the sense of (II)) of the system M; this abstract 
block diagram is then placed in the canonical form described in 
Theorem (4) of (II). We then adjoin to this abstract block diagram 


a family of mappings ®,, such that the following conditions are 
satisfied: 
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M.R.1I: The correspondence f —> ®, between the mappings in 
the abstract block diagram of M and the members of the 
family {®,} is 1-1 onto. 

M.R. II: If, in the block diagram of M, we have feH (A, A, ), 
then the associated mapping ®, is an element of the set 
n( TT B; H(A, A; ) , where the sets B; represent 

J 
the appropriate sets in the family © of environmental 
outputs of M. 

M.R. Ill: If ®, is the mapping corresponding to f, then fer (®;). 
The first of these conditions merely states that there is a mapping 
®, for every mapping f in the abstract block diagram of the system 
M. The second and third conditions assure that ®, actually be- 
haves so as to produce a replica of f from the appropriate outputs 
of M. 

We now proceed to discuss some of the general implications of 
the formalization developed above. 

A mapping f such that r(f) consists of a single element will be 
called a constant map. It is not hard to see that an arbitrary sys- 
tem M, all the mappings of which are constant maps, can display 
no fluctuation of outputs in a changing environment; i.e. are com- 
pletely invariant in their behavior. Therefore, the more closely 
the mappings of a system M approximate constant maps, the more 
invariant will its behavior appear over a wide range of environ- 
ments. It seems to be precisely the non-constancy of many of the 
mappings which occur in biological systems which accounts for 
their characteristic flexibility in dealing with environmental changes. 
Of particular interest in this connection is the possible non-constancy 
of the mappings which we have denoted by ®,, and which may be 
considered as representing a portion of the ‘‘genetic’’ material of 
the abstract biological system under consideration. A non-constant 
mapping ®, may under certain circumstances produce as output a 
mapping different from the mapping f with which it is associated, 
thereby perhaps altering the structure of the entire system. The 
implications of this type of situation, which is of obvious interest 
and importance, will be more thoroughly explored in another place. 

Thus far, our discussions of (M, R)-systems have not included 
any mention of the possibility of the replication of the components 
we have labeled FR, in (I), and denoted by mappings ©, above (cf. 
however, Rashevsky, 1958, in this connection). Nevertheless, such 
replication is one of the fundamental characteristics of living sys- 


116 ROBERT ROSEN 


tems, and must certainly be incorporated, in not too violently ad 
hoc a fashion, into any theory which makes a claim to general bio- 
logical significance. The representation theory which we have out- 
lined above enables us to make a number of remarks which may be 
pertinent to this problem. 

Quite generally, let X, Y be arbitrary sets, and let H(X, Y) de- 
note (for the moment) the full set of mappings f:X — Y. If we 
now fix an element zeX, and allow f to vary through H(X, Y), we 
find that we have an induced mapping 


b HIX, ¥y se 


defined by writing w, (f) = f(z) for all feH(X, Y). In this manner, 
we may make correspond to each weX a mapping yw, eH[H(X, Y), YI. 
It is immediately verified that for any two elements i,t, €X, Ve, = 
Ww, if and only if «,=2,; thus in effect we have constructed 
an embedding of the arbitrary set X into the set of mappings 
H{H(X, Y), Y]. It may be noticed that this procedure is a generali- 
zation of the standard embedding of a vector space in its second 
conjugate (or dual) space. 

Once again quite generally, let S, T be arbitrary sets, and let 
f:S —> T be a mapping in H(S, T). The mapping f then induces an 
equivalence relation fk, in S, defined by writing a, R,e, in S if and 
only if f(,)=f(#,). If we define the quotient space S/R, of § 
under 7 to be the set of equivalence classes of S obtained in this 
manner, we can form the diagram of mappings 


S—> s/R, 


Nh 


T 


where f is the given mapping, 7z is the natural mapping which maps 
each element of § on its equivalence class in S/R,, and f is the 
(uniquely determined) mapping which makes the diagram commuta- 
tive, It follows from our construction that fis a 1-1 mapping, and 
hence it has an inverse, f~!:T —» S/R, In particular, if f was 
originally 1-1, then S/R,;=S8 and f=f. Likewise, if the equiva- 
lence class of an element z,¢S consists of w, alone, then we may 
write f~' (f(#,)) =a, inS itself, 
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Let us now apply these general considerations to the abstract 
block diagrams we have defined above. Suppose that as before we 
enumerate the environmental outputs of the abstract block diagram 
of a system M as B,, B,,...,B,. As we pointed out above, the 
domains of the mappings ®, consist of the various cartesian prod- 
ucts ITB, . Suppose that f: A —»B is a mapping in the abstract 
block diagram of M, and that JJ B,, is the domain of the associ- 

Fé 
ated mapping ®,. If we now a X = ITB, ., and Y = H(A, B), 
then from the above general Sretectiuns we obtain an embed- 
ding of I B;. into the set of mappings 
je 


H = HHT ITB, H(A, 3) H(A, a} 


jel 


Now writing S = u( I By» H(A, B)), T = H(A, B), we see that for 
jeJ es 
each mapping Fe(H), we obtain an induced mapping F~', which 
maps the elements of H(A, B) into equivalence classes of map- 
pings in the set H/ J] Bi ys H(A, B)/e In case the equivalence 
jel 
class of a particular mapping consists of a single element, we have 
observed that the image of such a mapping under F~' may be con- 
sidered as an element of H/ I] Bia H(A, B)\. 
jes 
Let us now recall the significance of the various constructions 
we have performed above. Our embedding of the set TT B,, into 
jel 
the set of mappings (1) shows that a particular family of outputs of 
a system M/ i.e. an element in the set IT B; 


jeJ :) 


self be considered as a mapping. The domain of this mapping is 
the set a( HT B; H(A, B)\, which consists of mappings of the 
jeJ 


may sometimes it- 


type we have designated as ®,; the range of this mapping is the set 
H(A, B), (which in particular contains the mapping f:4 — B 
which represents a component of M). This mapping, called F in the 
above discussion, then induces a mapping F-1, which maps each 
element of H(A, B) onto an equivalence class of mappings of the 
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form ®,. The above discussion shows when this class of mappings 
is reduced to a single element. 

We have shown that, in the presence of suitable outputs of the 
system M, the mappings f of an abstract (M,R)-system may act in 
such a manner as to generate their associated mappings ®,. Thus 
we see that a means of replication of the components Ff, of an 
(M, R)-system, very similar in nature to the one postulated sepa- 
rately by N. Rashevsky (1958) is already contained within the 
formalism we have developed above, and requires the introduction 
of no new assumptions. It should be emphasized, however, that 
the above discussion of replication was introduced primarily to 
demonstrate the scope of our abstract formulation, and not neces- 
sarily as a definitive presentation of the ideas outlined therein. 


Il. The Time Lags. We recall from our earlier discussions in (I) 
and (II) that every system inherently possesses two different types 
of time lags; namely, the operation lags, caused by the actual func- 
tioning of the components, and the transport lags, which arise from 
the delays occasioned by the moving of materials or other stimuli 
from the neighborhood of the components which produce them to the 
neighborhood of the components to which they serve as inputs. We 
have already noticed in (I) that both types of lag depend on both 
the structure of the given system and the nature of the environment 
in which the system is operating. The formalization we have de 
veloped above makes it possible for us to define these lags (in an 
(M,R)-system) in a more precise fashion, and facilitates our in- 
vestigation of the possible relation of the time lag structure to the 
behavior of actual biological systems. 

We shall be concerned below with abstract (M, R)-systems satis- 
fying the conditions laid down inSection above. Letus suppose that 
f represents a mapping in a certain abstract (M, R)-system; we have 
observed that the time lag of the component represented by f de- 
pends on the particular input on which f is acting. That is, each 
object in the domain of f (i.e. each possible input to the component 
which f represents) will be associated with some positive real num- 
ber, which corresponds to the operation lag of the component in 
question when the given object is the input to that component, In 
more formal terminology, we see that the time lag associated with 
the operation of a component is a mapping, the domain of which 
coincides with the domain of the mapping f which represents the 
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component, and the range of which is contained in the set R* of 
non-negative real numbers. If we denote this mapping by 1,, we 
have 


t,:d(f) > per 


Naturally, similar considerations apply to the operation lags of the 
systems Ff, and represented by the mappings ®, We shall allow 
the value ~ as an admissible operation lag; as mentioned in Sec- 
tion I, a useful artifice for formally extending the domain of a map- 
ping f to a larger set A D d(f) is to replace f by a mapping f/f’ such 
that d(f’) = A, with the property that f(z) = f’(x) for zed(f), and to 
write r,(x) = r,-(x) for ved(f), ry, (@) = ~ otherwise. 

The transport lags are somewhat more complicated to represent, 
owing to the fact that, speaking roughly, they depend in general 
both on the component of origin and the component of destination, 
as well as on the particular material or stimulus involved. The 
most natural way of dealing with this problem seems to be the fol- 
lowing: let F be the cartesian product of the set of mappings of the 
system M of an abstract (M, R)-system with itself. For each inde- 
composable set A«M, we consider the cartesian product F x A (i.e. 
the set of all triples of the form (f, g, 2), where f, geM, aeA). We 
define the transport lag eve of the elements of A between the com- 
ponents of M represented by the mappings f, g respectively to be a 


mapping 


op te Rcd wee fad 


subject to the following conditions: 

1. aft =o if A is not a factor of d(g). 

2. a,/€ =o if A is not contained in r(f). 

3. o,f = 0 if (1) or (2) does not hold. 
The first two conditions express the fact that if the elements of the 
set AeM fail to be either inputs of the component represented by g 
or outputs of the component represented by /, then the transport lag 
between f and g of the elements of A is not defined (i.e. is infinite, 
for formal purposes); the third condition expresses the fact that the 
transport lag of the output of a component to the component itself 
is zero. With these conventions, we obtain a general definition of 
the transport lag structure of abstract systems which is consonant 
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with physical intuition, and which will be flexible enough for later 
applications. 

Let us now attempt to investigate in what manner the study of 
the time lag structure of an abstract (M, R)-system may throw light 
on other aspects of the system’s behavior. In particular, we desire 
to extend our discussion of the relations between the time lag 
structure and the ‘‘morphology’’ of abstract systems which was un- 
dertaken in (1). In order to do this, we must first make some as- 
sumption which serves to connect the two aspects of the structure 
of these systems. There is, naturally, a wide variety of possible 
assumptions which can be made in this context, each one corre- 
sponding to what N. Rashevsky (1956) has called a ‘‘system of ab- 
stract biology.’’ Hence it should be recognized that the tentative 
assumptions we put forward below are meant primarily to illustrate 
the type of results which may be obtained from these lines of thought, 
and therefore should not be considered as final. 

We make the following assumptions: 

L1: The transport lag oft is directly proportional to the physi- 
cal distance between the components represented by the 
mappings f and g; the constant of proportionality may vary 
with the particular ae A involved. 

In symbols, if we write p(f, g) for the distance between the com- 
ponents represented by f and g, this hypothesis becomes 


o,f® (a) =k(a)-p(f, 9). 


Of course, this hypothesis is meaningful only if the transport lag 
ot is defined (i.e. finite). 

I Psa i oie and o fe 8a are both finite, then the elements of A 
(which are outputs of the component represented by f) always 
go to the component, the representative mapping of which 
(either 7, or g,) gives the smaller of o i831, o fa . 

Physically, the assumption (L2) may be expressed as follows: the 
finiteness of o,f and o,f means, by definition, that the set A 
serves as a factor of both d(g,) and d(g,). By Assumption (Z1), 
these lags are measures of the distances p(T, 91). p(f, Gq), since 
in this case the same output a is being delivered by f to both 91 
and g,. Therefore, Assumption (L2) in its widest form requires 
that the output of a component of an abstract system is delivered in 
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full to the nearest component which is capable of accepting it as 
an input. 

L3: (‘Principle of Optimal Design’’) If f, 9,, g,,..., 9, are 
mappings of an abstract system such that the sets r(g,) are 
contained in the factors of d(f), and each factor of d(f) con- 
tains one of the r(g;), then the components represented by 
these mappings are spatially distributed in such a fashion 
that 


Bel bash 


— me es o ent 
r(&,) (8) . 


= rE) 

This principle specifies that if the mappings of a system are non- 
contractible (as we always assume; see (I)), then there are no ex- 
traneous ‘‘waiting lags’’ caused by the failure of inputs to arrive 
‘fon time.”’ 

Assumption (L3) is a precise formulation of the notion of op 
timality which was discussed in a more intuitive fashion in (I). 
Moreover, it seems to be a natural extension to cellular systems of 
the ideas of optimality developed by N. Rashevsky (‘‘Principle of 
Maximum Simplicity’’) in his discussions of morphology (1948), and 
considerably extended by D. L. Cohn (1954, 1955). 

The assumptions (L2) and (L3) are no doubt excessively restric- 
tive, but we remind the reader that they have been introduced pri- 
marily for illustrative purposes. Thus, for example, a more realis- 
tic assumption along the lines of (L2) would be to assume that, if 
a component represented by a mapping f has the property that r(/f) 
is a factor of the domains d(g,) and d(g,), where g,, g, represent 
other components of the system, then there is a distribution of the 
output of f proportionately between g, and g, depending on their 
relative distances from f. Likewise, a less stringent condition of 
type (L3) would require, instead of equality between the transport 
lags O, , that arbitrary differences of these lags should be 
smaller in absolute value than some preassigned positive number, 
which depends in general on the mapping f,; i.e. 


1? 


| Biof 


Ejoh 
ats aed eet) 


7(8;) 


for each i and j. For the present discussion, however, we shall re- 
tain (L2) and (L3) in their present form. 
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Let us now investigate in what manner the assumptions (L1)= 
(L3) affect the structure of an abstract (M,R)-system. Let us 
therefore suppose that a definite such system, which we may de- 
note by A, has been given. Then according to our previous discus- 
sions, we can find an abstract block diagram for A, which we shall 
denote by the same symbol A, and which consists of a certain num- 
ber n of mappings of type f, and hence also n associated mappings 
of type ®,. Since the range of each mapping (of either type) is an 
indecomposable set of the category (see (Il)) in which our repre- 
sentation of A is carried out, there must be at least n indecom- 
posable sets in the category (excluding the indecomposable sets 
which represent the purely environmental inputs to the system, and 
which are assumed to have zero transport lag). 

We observe that if there exists an indecomposable set A; 2 r(f) 
(for some fe A), such that A; . is a factor of the domain of two dis- 
tinct mappings g,, g,¢A, then (L1) and (L2) imply that, for the 
system to function properly, we must have 


e(f, 9,)=p(f, 9a) 


Otherwise, the total output of f will go to the ‘‘nearer’’ of I1> Go; 
thus the ‘‘farther’’ of g,, g, will by non-contractibility produce no 
output, and ultimately the entire dependent set of this mapping 
(see (I)) will fail to be produced. 

As a simple example of the type of conclusions which may be 
drawn from this last argument, let us suppose that A, represents 
an environmental output of M such that A; is a factor of the do- 
main of each mapping ®,, feM. Suppose that A; . > r(f), say. Then 
by the above, we must have 


e(f, ® ,) = constant 


for each feM. Recalling the physical significance of the objects 
involved, this condition signifies, roughly speaking, that the com- 
ponents corresponding to the mappings ®,(i.e. to the ‘“‘synthetic’’ 
or “‘repair’’ structure of the system A) is spatially distributed in a 
spherical fashion about the component represented by the mapping 
7. If we were to attempt to regard a free-living single cell as be- 
ing, to a first approximation, a system of (M, R)-type, then it would 
be natural to associate the mappings feM with biological functions 
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carried out by what is normally called the cytoplasm, while the 
mappings ®, would be associated with nuclear biological func- 
tions. If the cellular system also obeyed the laws (L1)-(L3), then 
the structural implication for the cell would be that the nuclear ma- 
terial (as observed cytologically) is distributed spherically around 
what would ordinarily be regarded as a purely cytoplasmic com- 
ponent. Hence this model implies an intimate relation between 
what we habitually call ‘‘nucleus’’ and ‘‘cytoplasm’’; in terms of 
relative biological activity, our model suggests the presence of 
certain biological activities which would usually be associated 
with cytoplasm, lying within the cytologically observed region as- 
signed to the nucleus. It is in fact well known (cf. for example 
Allfrey, Mirsky & Osawa, 1957) that isolated nuclei do carry out 
such cytoplasmic activities as the synthesis of proteins and the 
generation of high-energy phosphate bonds. 

Let us now carry the above construction one step further. Sup- 
pose as above that A; represents an environmental output of M 
which factors the domain of each mapping ®, of the system A, 
Suppose in addition that there exist environmental outputs 4;,, 
A;,,..., Ai, of M such that, for each mapping ®;eA, exactly one 
of the sets A;,, Aig,---, Ai, factors d(®,). This process in ef- 
fect decomposes the set of mappings [Drbrem into r classes, the 
mappings in the z*® class being such that A;, is a factor of their 
domains. From what we have already said above, it follows that 
the mappings in the &* class must lie on a sphere S,, such that 
the component which produces the output of M represented by the 
set Ai, lies at the center of S,. But we have seen that every 
mapping ®,, regardless of its class relative to the sets Ai,, Aig, 
..., Ai,, must already lie on the sphere S, the center of which is 
the component which produces the output A;, (since by hypothesis 
the set A;, factors the domain of every mapping ®,). Hence the 
mappings ®, in the 7*+ class are restricted to lie on the intersec- 
tion of the sphere S, with the sphere S,- This intersection 1s a 
circle; the important thing to notice, however, is that the circle is 
a one-dimensional manifold. The fact that this model requires the 
‘“genetic’’? material of the system A to lie along one-dimensional 
manifolds is suggestive of the linear arrangement of the hereditary 
material found uniformly in all types of organisms. Our model al- 
ready incorporates the feature that the mappings ®, are not ran- 
domly distributed along the sphere S,, but are organized linearly 
according to their inputs and outputs according to the rules (L1)- 
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(L3); this structure may be related to a conjecture of M. Demerec 
and Z. E. Demerec (1956). These authors suggest that the linear 
arrangement of the genetic material of cells may be connected with 
the existence of a sequence of coupled chemical reactions, each 
of which is mediated by an enzyme controlled in turn by a particu- 
lar gene in the linear array. The reaction sequence is arranged in 
such a manner that the 7*® reaction of the sequence is mediated by 
the enzyme controlled by the i‘® gene in the linear array. We shall 
not, however, pursue this point any farther at this time. 

It should finally be noticed that the basic characteristics of the 
above model are not greatly affected by a weakening of conditions 
(L2) and (L3) along the lines mentioned previously on p. 121. The 
reader may readily convince himself that the principal effect of 
such a weakening will be to introduce a ‘‘spread’’ about the inter- 
sections of the various spheres which we have designated by 83,5 
S;, above; within reasonable limits, this spreading of circles into 
tori does not affect the conclusions we have drawn from the model. 

These simple results already indicate the directions in which a 
theory of the type developed above may lead. As remarked in (1), 
the segregation of the synthetic material into nuclei by living cells 
is one of the most commonplace observations of biology, and yet is 
one for which no theoretical justification of any kind has been pro- 
vided. We emphasize once more that the above very simple con- 
siderations should not be taken too seriously in this regard; they 
were presented primarily for their heuristic value, in pointing out 
the possible usefulness of some of the methods developed above. 
However, the fact that our results are for the most part invariant 
under a loosening of some of our hypotheses may indicate that by 
imposing more realistic conditions on the time-lag structure of 
(M, R)-systems, we may obtain results in closer accord with ob- 
servation. 

Let us now turn our attention to certain other aspects of the 
physical structure of (M,R)-systems, which may be inferred from 
the hypotheses (L1)-(L3). We propose to investigate the relations 
between the abstract block diagram of the (M,R)-system A and the 
number of its possible different optimal realizations. To begin 
with, let us consider a single mapping of type DeA and let us 
suppose that 


d(®,)=A,;, x A;, x 10 X Aj, 
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and Aj, s e(fi,)- tor k= 1, >... $8; where the fi, are terminal map- 
pings of M. If the entire Sere A begins to operate at time ¢ = 0, 
we shall denote the time at which the first element of the set Ai, 
produced by the system appears as T;, (thus T;, is, roughly SO6ale 
ing, the time elapsed between environment and output). Each 7;, 
is readily seen to be a sum of operation and transport lags in M. 
The non-contractibility of the mapping ®, implies that the com- 
ponent which corresponds to ®, cannot begin to operate (to produce 
a replica of the component corresponding to the mapping f, or some 
related structure) until the time 

iy fi,» % 


. 1 f fi ’ 
max (Ti, + oA; ) Pigt oA, gl eeey Ce 


If we now assume that the system A is optimal (i.e. that the condi- 
tion (L3) is satisfied), we must have 


=...=7, +0,/'8 
ty lg s Us 


This is a set of s-—1 conditions on the operation and transport 
lags of A. 

Thus we see that, using the hypothesis (L3), the abstract struc- 
ture of the block diagram of an (M,R)-system, we impose condi- 
tions upon the time lag structure of the system. Let 9 represent 
an arbitrary mapping of an (M,R)-system (i.e. of type f or of type 
®,). We shall denote by g(?) the number of factors of the domain 
of ~, which do not represent purely environmental inputs (i.e. fac- 
tors which contain the ranges of other mappings of the system; 
these are the only types of set which impose conditions on the time 
lag structure). Then we observe immediately that precisely the 
same argument as given above implies that for each mapping feM, 
the mapping ®, implicitly imposes a total of g(®,)-1 conditions on 


A, and hence there are Ys (7(®,) -1)= Di q(P,)-n conditions 
fem fem 

on the lags of A, due to the mappings @, alone, (where n is the 

number of mappings in M). To enumerate the number of conditions 

placed upon A by the mappings feM, we observe that for the map- 

pings /, as for the mappings ®,, 9(f) non-environmental factors in 
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d(f) implies a total of g(f) — 1 conditions on the lags of M. Hence 
the assumption of hypothesis (L3) imposes a total of 


Y. G@)-+ Di @A-H= D, ap) +9) - 20 


fem fem fem 


conditions. 

In a somewhat similar fashion, the hypothesis (L2) also imposes 
a number of further conditions on A. More precisely, let A be a set 
in A which represents a non-environmental input to some component 
(i.e. A> r(y) for some mapping weA). Let us denote by p(A) the 
number of mappings (eA such that A is a factor of d(p~). Then the 
same line of reasoning which we have employed above shows that 
there must be a total of p(A) — 1 conditions imposed on A by each 
such set A, and hence there will be 


>. (P(A) -1) 


AeA 


conditions on the lags of A by the hypothesis (L2). Since there 
must be n such sets, we find that this number can be written 


a 
2d, P(A)=n. 


AeA 


Finally, there are n further conditions to be imposed on the lags of 
A if we make our usual assumption concerning the finite life-time 
of the components of A (see (1) for the details), Combining all of 
these conditions, we find that we have a total of 


2G [a (f) + 9(®,)] + >; p(A)-2n 


feA AeA 


separate conditions on the time lags of A. 

It is clear that if n is the number of mappings in the system M of 
A, then there are a total of 3n lags in A, of which n are transport 
lags, and 2n are operation lags due to the action of the compo- 
nents represented by the mappings f and ®,. Combining this in- 
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formation with what we have previously obtained, we find that we 
are confronted with two possibilities: 


tte >: [a(f) + g(®)] + ms p(A)-2n<38n; 
feA AeA 


2. D) (g(f)+9(@)l+ )” p(A)-2n>3n. 
fed AeA 


If the first possibility is realized, then we have too few conditions 
to completely determine the lag structure, and hence there are 
many (in fact, infinitely many) different optimal forms which the ab- 
stract (M,R)-system A can assume. Furthermore, since the given 
conditions only tell us how many of the lags we may solve for, but 
not which ones, it is entirely arbitrary which lags we choose to 
solve for, and which ones we choose to select at random. Hence 
we discover a great latitude in constructing different optimal forms 
of the same abstract system; we may choose to completely fix the 
operation lags of the system (if we have a sufficient number of 
conditions), thereby obtaining systems which differ in the relative 
distances between components, or we may choose to leave certain 
operation lags undetermined, thereby opening up the possibility of 
choosing different physical realizations for the same component 
(cf. our discussion of ‘‘amplifier’’ in (1), p. 247). Thus, a study of 
what we may call the ‘‘morphology”’ of abstract (M, R)-systems of 
this type is of a degree of complication approaching what we have 
become accustomed to in the study of actual biological systems. 

If the second possibility holds, we face the existence of two 
further sub-possibilities; namely, either all the conditions imposed 
on the system are all independent, or else they are not indepen- 
dent. If the conditions are independent and the equality sign holds, 
then there is exactly one optimal form for the system in question. 
If the conditions are independent and the inequality holds, then the 
system in question cannot be put into an optimal form. If the con- 
ditions are not independent, then the total number of independent 
conditions among them must total less than 37 in order for optimal 
systems to exist in great numbers. In systems of the type dis- 
cussed above, it is immaterial whether the imposed conditions are 
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independent or not, since the totality of conditions is already less 
than the number of lags of the system. 

Thus we see how the principles which we have advanced above 
may enable us to make inferences concerning the structure of 
physical realizations of abstract systems, using the abstract struc- 
ture alone. As we have emphasized above, it is not the specific 
form of the principles which we have labeled (L1)-~(L3) above 
which is of primary interest; it is the fact that such principles can 
lead to useful results that is of significance. We emphasize again 
that the above discussion is presented primarily as an example, to 
illustrate in what manner similar, but more realistic, principles can 
be used to obtain results of a theoretical nature which may have 
general biological importance. 


The author is indebted to Professor N. Rashevsky for a thorough 
discussion of the manuscript. 

This work was aided by United States Public Health Service 
Grant RG-5181. 
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Diffusion through a flat pore into a large open region is proportional 
to the linear dimension of the pore and not to its area. This was first 
explained by Brown and Escombe (1900) for a circular pore and is here 
generalized, by means of a dimensional argument, to include any type of 
regular opening. The problem is further generalized to include diffusion 
through pores of finite thickness, finite distance apart, and into finite 
regions. Since this problem cannot be solved exactly, an approximation 
method is introduced. Reasons for the credibility of the approximation 
are presented. It is then shown, by means of the approximation method, 
that the diffusive flow through a pore is equal to the total concentration 
difference divided by the resistance of the system. The resistance, in 
turn, is the sum of the resistances of all portions of the system, each of 
which is calculated. The result is compared with results which have 
been calculated exactly for limiting cases and found to agree very well. 
The results are then applied to a standard method of computing pore 
size in membranes, and it is shown that the correction factor is negli- 
gible. 


Although the open surface area of the stomata of plant leaves is 
of the order of 1 per cent of the total leaf area, the evaporation of 
water from the leaf is of the order of 50 per cent of the evaporation 
of an exposed surface, the area of which is equal to that of the 
leaf area (Meyer and Anderson, 1952). This was explained by 
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Brown and Escombe (1900), who showed, both experimentally and 
theoretically, that the diffusion rate of a neutral substance through 
a flat hole (i.e., a pore with zero thickness) into an infinite vol- 
ume was proportional to the diameter of the hole and not to its 
cross-sectional area. Thus, for the leaf, where the stomata are 
sufficiently far apart that the presence of the other stomata inter- 
feres little, if any, with the flow through a given stoma; the flow 
through one large hole is less than the total flow through many 
small holes, the total area of which is equal to that of the large 
hole. This is because the diameter of the large hole is less than 
the sum of the diameters of the small holes and the flow is pro- 
portional to the diameter. 

For the flow through a hole which is regular in shape and can 
be characterized by a single variable (such as the radius for the 
circle), the dependence of this flow on the variable may be de- 
rived by means of a simple dimensional argument. Let L represent 
length, T time, and M mass. If the hole is of zero thickness, if 
the diffusion takes place into a uniform infinite medium, and if the 
material which is diffusing has a concentration C, at the hole and 
a zero concentration at infinity, then the only variables are j, the 
total flow through the hole, with dimension MT™!; D, the diffusion 
coefficient, with dimension L?77}; C,, the initial concentration, 
with dimension ML~*; and a, the linear dimension of the hole, with 
dimension L. Denote by & a dimensionless constant; then the only 
dimensionally consistent relation between these four quantities is 


j= kDCya. (1) 


For the case of a circular pore, equation (1) is in agreement with 
the results of Brown and Escombe (loc. cit.) with & equal to 4 and 
a equal to the radius of the pore. 

Almost all diffusion problems found in the literature deal with a 
flow proportional to the area through which the material is dif- 
fusing. Therefore, it may be asked why in this case the flow is 
proportional to the diameter of the pore. The answer appears to 
be that the major resistance to the flow is offered not by the pore 
_ (since the latter is of zero thickness and thus has zero resistance) 
but by the infinite region above the pore. Thus, a priori, the flow 
could be proportional to any power of the diameter. When a ma- 
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terial passes through the center of a hole, symmetry conditions de- 
mand that it can flow only in an axial direction, whereas when it 
passes near the perimeter of the hole, it can flow axially as well 
as in a direction parallel to the hole. Therefore, the effective re- 
sistance which the infinite region offers to different parts of the 
hole differs, the center having a higher resistance than the edges. 
Thus the flow along the sides of the hole will be larger than the 
flow through the center. For Poiseuille’s law (flow proportional to 
a*) the flow along the edge is zero. For the usual diffusion law 
(flow proportional to a?) the flow along the edge is the same as 
the flow through the center. Thus it is not unreasonable to expect 
that the flow through a flat pore into an infinite region, where the 
flow along the edges is greater than the flow through the center, 
should be proportional to a. 

The above discussion is applicable only to a very limited situ- 
ation, namely, a pore of zero thickness and infinite regions around 
it. The case of a pore of finite thickness but infinite regions 
around it has been treated by Brown and Escombe. The case of 
finite regions and pores of zero thickness has been treated by 
Verduin (1949) by a different approach from that used here. A more 
general picture in which there are a number of pores of finite 
thickness, with a finite region on either side of the membrane, will 
be discussed in this paper. The usual simplifying assumptions 
will be made, namely, that the pores are of uniform size and thick- 
ness and are uniformly distributed within the membrane. 

To the best of the author’s knowledge, this case cannot be 
solved in general with present-day techniques. It can be solved 
only for specific cases by numerical methods. General solutions 
can be found only for the case discussed above (a zero-thick pore 
with infinite region around it) or for the case where the radius of a 
zero-thick pore is very much smaller than the height of the cylin- 
drical region above the pore, which height in turn is very much 
smaller than the radius of the cylindrical region (Gray, Mathews, 
and MacRobert, 1931). For this reason, an approximation approach 
should be tried. The particular method used here is similar to 
that first proposed by Landahl (1953). 

In the approximation method used in this paper, we assume that © 
circular pores of thickness w and radius a are uniformly dis- 
tributed within a membrane, with an average distance 2d between 
the centers of adjacent pores (see Fig. 1). Thus each pore has a 
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FIGURE 1 


cylindrical region of radius d which ‘‘belongs’’ to that pore. Ata 
distance S, from the left side of the membrane there is a mixed 
region in which the concentration of the diffusing material has the 
fixed value C,; at a distance S_ from the right side of the mem- 
brane there is another mixed region in which the concentration of 
the diffusing material has the fixed value 0. These assumptions 
are not too drastic and are in line with previous attempts to treat 
this problem. However, the following assumptions will be much 
more severe. 

It is assumed that the flow from a pore passes first through a 
funnel-shaped region (see discussion below for rationale). When 
the larger opening of the funnel reaches the point where its radius 
is d, it is assumed that the flow is then confined to a cylinder of 
radius d. For the flow into a pore (from the left side) the se- 
quence is reversed. (The funnel and cylinder are shown as dashed 
lines in Fig. 1.) The flow through the pore itself will be as if 
through a cylinder of radius a. The other major assumption is that 
the concentration of diffusing material is uniform for any cross- 
section parallel to the membrane. The rationale behind these as- 
sumptions will now be discussed. 

The assumption that the flow from the pore is initially through 
a funnel is a compromise between the assumption that there is no 
spreading-out from the sides of a pore (the flow is through a cylin- 
der of radius a) and the assumption that there is complete spread- 
ing-out from the sides of the pore (the flow is through a cylinder 
of radius d). Also, for the approximation to have any plausibility, 
the slope of the sides should be of the order of unity. (Later in 
this paper, the slope will be chosen to be 4/z, which is close to 
unity, in order to have exact agreement between the results of this 
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approximation and the true results for the case of S, and w equal 
to zero and S_ and d equal to infinity.) 

The resemblance to a funnel close to the membrane and to a 
cylinder for larger distances from the membrane is also obvious if 
one examines the flow lines for the above case, as drawn in Brown 
and Escombe (Fig. 7 of their paper). The limitation that the dif- 
fusion must first occur through the funnel region restricts the total 
amount of space available for diffusion and therefore decreases 
the diffusive flow. However, the assumption of a uniform concen- 
tration for any cross-section parallel to the membrane implies that 
the diffusion coefficient for diffusion in a direction parallel to the 
membrane is infinite. This tends to increase the flow over that 
which would be present through the funnel without this assump- 
tion. Therefore, the two assumptions affect the total flow rate in 
opposite directions, and by a suitable choice of the slope of the 
funnel sides, a good approximation to the true flow may be 
reached. 

The major criteria of the goodness of this approximation method 
is, of course, how well it agrees with the exact result. The only 
cases which can be solved exactly, besides the two previously 
mentioned, are for both S’s to be infinite (in which case the flow 
would approach the flow in a cylinder of radius d) or for S to be 
zero (in which case the flow approaches that through a cylinder of 
the size of the pore). As seen from the assumptions listed previ- 
ously, the approximation method for these two cases yields the 
same results. Finally, for the previously mentioned case in which 
a is much less than §, which in turn is much less than d, and 
where S) and w are zero, the exact results turn out to be the value 
for the case of the infinite region plus correction terms. The ap- 
proximation method yields the same type of result with the cor- 
rection term agreeing to within 14 per cent error. Thus, on the 
basis of all the cases for which there is an exact solution, the 
agreement is as good as could be expected. Also, as discussed 
in Landahl’s paper (loc. cit.), this type of approximation approach 
is in good agreement with the exact solutions for almost all the 
problems to which it has been applied. 

An exact solution would obviously not apply in the following 
situations: where the pores are not perfect circles; where there is 
not a uniform distribution of pores, so that the assumption of a 
cylinder of radius d about each pore is not valid; and where the 
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pores do not go straight across the membrane but may be slanted 
and not of uniform cross-section. However, the approximate solu- 
tion should still apply as well as before. (I.e., a slight change in 
the geometry of the problem will drastically alter the exact method 
of solution as well as the form of solution, but it probably will 
not appreciably affect the general form of the approximate solution 
and the values of its parameters.) 

Finally, it should be pointed out that, although Brown and 
Escombe stated that they were solving their problem exactly, they 
actually used, in calculating the flow within the pore, the assump- 
tion of a uniform concentration for any cross-section parallel to 
the membrane. 

It should be emphasized that the approximations used in this 
paper are at variance with those used by Verduin (Joc. cé¢.) in his 
attempt to solve a similar, but more restricted, problem. However, 
it should be noted that he does not take into account all the pa- 
rameters that should be considered (i.e., he ignores the S’s, w, 
and a) but incorporates them in his equation by means of a con- 
stant of proportionality and an arbitrary constant (which is evalu- 
ated for d=), Verduin also assumes that, since the diffusion 
from a single pore forms diffusion spheres, the concentration 
varies as the inverse square of the distance from the pore. In 
reality, it should vary only as the simple inverse of this distance. 
Further, he assumes that this distribution will still hold true when 
there is interference between pores. Because of these factors, 
the approach of Verduin will not be used. However, since the ex- 
perimental data do not contain all the parameters that are needed 
in order fully to test the results of the approach of this paper, it 
cannot be definitely stated which approach is better. 

The major method used to solve the specific problem here is 
based on the assumption of continuity. In the steady state the 
total amount of material crossing any cross-section parallel to the 
membrane is the same for all cross-sections of the flow. A similar 
method can also be used to solve the problem of the transient 
state, but, because of the limited value of this solution and be- 
cause of the extreme algebraic complexities involved, only the 
steady state will be considered in this paper. 

Let « represent the distance. from the left-hand boundary of the 
system to any cross-section parallel to the membrane within the 
system; A(z) be the area of the cross-section parallel to the mem- 
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brane at the distance 2; j be the total flow through the pore per 
unit time, which in the steady state will be the total flow across 
any cross-section parallel to the membrane per unit time; C(z) be 
the concentration of the diffusing material at the distance x; and 
D(z) be the diffusion constant at the distance z. Then, owing to 
the assumption of uniform concentration for any cross-section 
parallel to the membrane, 


dC(z) 


j = — D(x) A(z) oe 


: (2) 


Since 7 is a constant, the above equation can be readily inte- 
grated. Denoting by / and r the left and right sides of the system, 
we obtain 


3 os 
j=-——2°_., (3) 


r dz 
J D(z)A(z) 


Using the analogy of an electric circuit, j corresponds to the cur- 


Tr 
rent, C,, to the E.M.F., and therefore | [1/D(x)A(2)ldz corre- 
l 


sponds to the resistance Ff. Thus, since the integral from 7 to ris 
the sum of the integrals over all the subregions into which the 
system may be subdivided, the total resistance of the system is 
the sum of the resistances of all these subregions. (It should be 
emphasized that the ability to break up the system in this manner 
and consider the total resistance to be the sum of the subresist- 
ances, with no interaction between the regions, is not valid for the 
exact case and arises solely from the assumption of a uniform 
concentration for any cross-section parallel to the membrane.) 

It is most convenient to subdivide the system into the obvious 
regions of the funnels, the pore itself, and the cylinders of radius 
d. Thus, if 1/m represents the slope of the sides of the funnel, 
and D(a) is constant within any subregion, it may be easily de- 
rived from equation (3), that the resistance, F,, for the cylindrical 
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region is 


= 5 a ean S >m(d - a), (4) 
s Dd? 
the resistance, Kies for the funnel is 
(d — 
gine he tn 1 esl pear (5) 
f D, nad 


while the resistance, H,, for the pore is 


Res (6) 


The results for R, and R, are what would be expected for the dif- 
fusion through a cylinder. For the case where S is not large 
enough for the radius of the funnel to reach the value d, the funnel 
resistance becomes 


S 
Pe aS aera. (7) 
D,nwa(am + 8) 


In order to calculate m, it is necessary to consider the case where 
d is larger than S, S tends to infinity, and w is zero. For this 
case, as stated earlier in the paper, the actual total flow through a 
pore is 4DC,a. According to equations (3) and (7), the flow is 
given by DC,na/m. Therefore, 


(8) 


a) 


As stated before, the only other exact case treated in the liter- 
ature (Gray, Mathews, and MacRobert, loc. cit.) is for w =0, and 
a<<S<<d, The exact solution yields a total resistance equal to 
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zat i Sea pS 
4Da aDS DS 


The approximate solution for this case, equations (7) and (8), 
leads to a resistance equal to S/[4aD(S + 7a/4)]. This expression 
may be expanded in a power series, and, since (a/S) is taken to be 
very small, only the first-order terms of (a/S) will be retained. 
Thus the resistance is 


A 7 
4Da 16DS° 


The first term is in agreement with the first term of the exact re- 
sult. Since (ln 2)/m7 = 0.222, while 7/16 = 0.195, the second terms 
differ by 14 per cent. The third term of the exact solution is, for 
d = 8, 1/56 of the second term and is much less for d larger than S. 
Thus the omission of the third term is not a bad approximation. 

The above results may now be applied to the study of diffusion 
of water through animal membranes. For most, if not all, of the 
cases reported in the literature, the flow through a membrane plus 
a thin region of water on either side of the membrane is much less 
than the flow through a correspondingly thick sheet of pure water. 
Since the flow through the sheet of water is the same as through a 
membrane system with the pore radius equal to d, d must be much 
larger than a. Moreover, the resistance to the flow of water 
through the membrane due to the region of the cylinders of radius d 
on either side of the membrane is negligible as compared to the 
resistance of the flow due to the funnels and to the pore itself. If 
the further reasonable assumptions are made that § is much larger 
than a and that D is the same for all the subregions, then, because 
of equations (3)-(8), the flow through a pore is 


2 
atl DrC ya (9) 


7 
wt+-a 
2 


This is the result of Brown and Escombe for the more limited case 
of diffusion through a pore with d>>S—>+~. Note that for this 
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case the major resistance of the regions outside the membrane is 
due to the small radius portion of the funnel, i.e., the portion of 
the funnel close to the membrane. The above can be applied to 
the studies of the value of a@ for various animal membranes. The 
method used by Koefoed-Johnsen and Ussing (1953) and Pappen- 
heimer (1953) to find the ‘‘equivalent’’ pore consists in determin- 
ing the flow of water through a membrane by means of diffusion 
and by means of an osmotic pressure difference. Then, assuming 
Poiseuille’s law to hold and taking the ratio of the two flows, the 
value of the radius of the pore may be calculated. However, they 
did not include the effect of the outside solution located in the 
funnel. .If this is included and if Poiseuille’s law is assumed, 
i.e., if the osmotic flow through a pore is proportional to a‘*/w, 
then it can be shown by simple algebraic manipulation that 


a=, (10) 


where a, is the value of a under the assumption that the regions on 
either side of the pore have no effect on the diffusion flow. Note 
that the ‘‘true’’ value of @ is always less than dy. Since it has 
been found that a, is less than 50A (more specifically, it is of 
the order of 5A) and since w is greater than 100A, the correction 
term is at most 30 per cent and is probably of the order of 5 per 
cent or less. Further, as pointed out by Dr. P. Curran (personal com- 
munication), Poiseuille’s law fails when a/w > 1, which is the only 
range in which the correction factor is appreciable. Thus the ef- 
fect of the surrounding region upon the diffusion flow through a 
pore in an animal membrane is not of much importance for the type 
of problems considered here. 

However, the experiment used to find the diffusion flow should 
be designed so that the resistance to this flow due to the sheets 
of water surrounding the membrane is small compared to the re- 
sistance due to the membrane plus funnel. For example, Koefed- 
Johnsen and Ussing (loc. cit.) found that the resistance to the 
diffusion of D,O across the frog skin is only four times the re- 
sistance of an equivalent thickness of water. Since the portion of 
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the skin which offers the major resistance to the flow of D ,J con- 
sists probably of one or two membranes, which are much Cine 
than the skin itself, the resistance of the skin is effectively the 
resistance of these membranes plus funnels plus the resistance of 
a sheet of water the thickness of the skin. The sum of these re- 
sistances being four times the resistance of the water sheet alone, 
the resistance of the membranes plus funnels is therefore only 
three times the resistance of the water sheet. Hence the diffusion 
flow across the membrane alone would be 4/3 of the diffusion flow 
across the whole skin. The same argument would also apply to 
the results of Durbin, Frank, and Solomon (1956) on frog gastric 
mucosa, with the factor 4 in the above case replaced by 3.6. Hence 
the diffusion flow across the membrane alone would be (3.6/2.6) of 
the diffusion flow across the whole mucosa. For this case, the 
equivalent pore radius would be decreased by a factor of 
/(2.6/3.6) = 0.85. 

For the diffusion of water or CO, through a plant stoma or 
through artifical membranes, the value of a/w may be of the order 
of, or greater than, 1, and the surrounding region is then of great 
importance. But, in order to test the above equations experi- 
mentally, all the parameters (S,, S,, w, and d) must be controlled. 
This has not been done in the experiments with which the author 
is familiar. 


The author wishes to thank Professors H. D. Landahl and Robert 
Macey and Mr. Stuart Zimmerman for their many helpful discussions 
and suggestions. 
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The transient stage of the random dispersal of logistic populations is 
investigated, using a Sturm-Liouville series leading to an infinite system 
of non-linear integral equations. These equations are then solved via a 
Successive approximation scheme. R. A. Fisher’s (steady-state) velocity 
of advance paradox is discussed. An illustrative example is worked to 
the second order of approximation. 


Introduction. The random dispersal of populations of living 
organisms subject to a logistic growth rate is a natural extension 
of the Lotka-Volterra theory to include space dependence, as well 
as time dependence, of the population density function. Unfortu- 
nately, the introduction of a logistic growth rate leads to a formida- 
ble non-linear partial differential equation describing the biologi- 
cal situation. There are two ways to avoid the difficulties of solv- 
ing the time-dependent equation. The first method merely neglects 
the time dependence of the density function by assuming that a 
steady state exists; the analogous genetic problem was first 
treated in this manner in a classic paper by R. A. Fisher (1937). 
In this paper, Fisher was led to a paradox concerning the inde- 
terminacy of the velocity of the ‘‘wave of advance of advantageous 
genes” because his analysis neglects the initial conditions. Simi- 
larly, J. G. Skellam (1951), treating the steady-state ecological 
problem, encountered the same difficulty. The steady-state proh- 
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lem is really anomalous because we are making a statement about 
the biological system that holds for all time. The Lotka-Volterra 
biological mechanics, however, is involved only with the predic- 
tive description of the system from specified initial conditions. 

The second method retains the time dependence but involves a 
perturbation scheme designed to linearize the equation by con- 
sidering a biological state near equilibrium. This scheme is also 
unsatisfactory, since it is valid for only very small time intervals. 

Finally, we are led to consider the equation without any simpli- 
fying assumptions. To the author’s knowledge, the only investiga- 
tion of the transient case has been by H. D. Landahl (1957), who 
used an interesting approximation technique. The transient-state 
analysis can also be attacked hy utilizing a Sturm-Liouville series 
leading to an infinite system of non-linear integral equations. In 
this manner, we treat the problem as an initial-value one and pass 
to the steady state hy letting the time become infinite. Using this 
technique, we can resolve the indeterminacy of the velocity by 
taking the initial conditions into account in the analysis. Although 
we can prove that there is some unique velocity, it is not neces- 
sarily equal to Fisher’s minimum velocity of the ‘‘wave of ad- 
vance.’’ The further advantage of this method lies in its construc- 
tional value; we can solve the resulting nonlinear system of inte- 
gral equations to any desired degree of accuracy in the space-time 
domain. 

Method of solution. In what follows, we restrict ourselves to the 
one-dimensional case; extension to the two-dimensional case is 
immediate, and we do not consider it in this paper. 

The usual one-dimensional equation describing the random dis- 
persal of a logistic problem is given by 


Oy =, oy 
—-a@ 
ot Ox? 


F y¥(1-¥y), (1) 


Where y(z,¢) = population density, a? = diffusion constant, and 
y = logistic growth coefficient. This equation is derived under the 
assumption that the dispersion is normal. If, however, non-normal 
dispersion were to take place, then our fundamental equation (1) 
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would take the form 


= 8 er ed 
ee av iy), (2) 


where K, is the nth cumulant of the non-normal distribution con- 
sidered. When normal dispersion occurs, then (2) degenerates 
to (1). 

Equation (1) is not the most general description of the normal 
process, for we can consider the diffusion and logistic coeffi- 
cients as being functions of distance instead of being merely con- 
stants. The basic equation then is 


on ) 


au 
=~ 5 [P(@)=-| = ye) ¥(1- wy) (3) 


Equation (3) is a non-linear partial differential equation of para- 
bolic type, the non-linearity being due to the presence of the func- 
tion y(z)¥w(1-wW). It is possible to solve equation (3) (under 
fairly mild restrictions) using a Fourier method (Sundaram, 1939), 
which consists of converting the partial differential equation into 
an infinite system of non-linear equations. This infinite equation 
system is then solved by using the method of successive approxi- 
mations; in this manner we account for both boundary and initial 
conditions explicitly. 

Let p(x) (distance-dependent diffusion function) be a positive 
bounded continuous function with bounded first and second de- 
rivatives: 


0S p(z) <M, lpta)l<m, le@)l<m. (4) 


We try to determine a solution of equation (3) which is regular in 
the domain 


Gs2% nO Se So, (5) 
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and which will satisfy the initial condition 
(2, 0) = (2) (0S 2 Sz) (6) 

and the boundary conditions 

w (0,t) = wy (7,t) = 0 (¢ 2 0). (7) 
(Without any great change in the argument we can consider other 
boundary conditions; we merely use the specific boundary condi- 

x ¥ 

tions to illustrate the method.) 


We write the solution of (3) as an infinite series of eigenfunc- 
tions 


U(a,t) ~ D° U,(t)}d, (2), (8) 


n=1 


where the ¢, (z)’s are the complete ortho-normal eigenfunctions of 
the Sturm-Liouville differential equation (inO<S zz), 


d dd 
with the boundary conditions, 


$(0) = d(m) = 0. (10) 


Denote by A,, the eigenvalue corresponding to d(x). We can prove 
that the eigenvalues A, are non-negative for the boundary condi- 
tions (10): 


aaa f [o(@]*ae~ [99 ae, (11) 
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where the first term vanishes and the integral is positive (except 
when ¢ is a constant), 
The initial-value function can be written 


o(z)~ >) C, $,(2), (12) 
n=1 


where C, = U_(0). Finally, we assume the expansion of 
f(2, t, vs) = y (x) (1 > vs) (18) 


in a series of ¢, (x). As Sundaram has remarked, a necessary con- 
dition that f(a, ¢, uw) satisfy the boundary conditions (7) is that 


(0, ¢, 0) = f(z, t, wv) =0. (14) 


Now let 7 have the expansion 


f(z, t, 4) ~ > F,() }$,(@)s 


n=1 


(15) 


F(t) = i ahold Spy ie: 


Our problem reduces to the determination of the coefficients U, (t) 
such that (3) and (6) can be satisfied. Substituting (8) in (3) and 
equating the Fourier coefficients ¢, (#) on hoth sides, we find that 
U_,(¢) must satisfy the non-linear differential equation system: 


dU, 
dt % An U,, = F,(), 


(16) 


U(0)=C, (n=1,2, +++) 
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(We could throw this into matrix form and consider the topology of 
the system, but this will not be done.) Equation (16) is equivalent 
to the non-linear integral equation system: 


t 
U. (t) = cet? ri [ go Aatt-s) F ls, U(s),...] ds. (17) 
0 


Sundaram (1937) has shown that this system can be solved uniquely 
by the method of successive approximations, provided that we im- 
pose the condition that 


yn On < (18) 
n=1 


in fact, it can be shown that the convergence of this series shows 
that w(x) exists and is integrable L*. These conditions upon the 
initial-value function w,(z) are not severe but do exclude certain 
functions which do not die out rapidly enough at infinity and yet 
are of possible ecological interest. 

If we use the successive approximation scheme, 


U(t) = 0, (19) 
t 
U“™ (t) = Ce fat / Pe hg FL fay UM (a) s, ol das 
0 
we can prove that 
; ) 
lim US™ (t) = U,, (t) (20) 


will exist for all ¢ in any finite interval and all m. 
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Velocity of the wave of advance. The solutions of equation (3) 
must be non-negative functions, seeing that they represent popula- 
tion densities. Under this restriction of positive-definite solution 
and assuming that a steady state exists, R. A. Fisher (1937) 
showed that progressive waves are possible with any velocity 
oreater than, or equal to, a certain minimum velocity determined by 
both the diffusion and the logistic growth coefficients where these 
coefficients are now taken as constants. (In the remainder of this 
paper we restrict ourselves to constant diffusion and logistic 
growth coefficients, not so much for mathematical reasons as for 
ecological reasons. We do not know enough about these coeffi- 
cients at present to warrant the added mathematical complexity of 
their space dependency.) 

If we restrict the logistic growth coefficient y to be greater than 
zero (growth problem), then we can show (Fisher, 1937; Kendall, 
1948) that the only possible velocities of propagation are 


lv|>2ay%. 


The minimum velocity ay’ will be termed v,. This velocity of 
propagation of the wave of advance is independent of the initial 
conditions and indeterminate because there is no one definite ve- 
locity of propagation but a multiplicity of possible velocities. If 
we take cognizance of the initial conditions, we should be able to 
show that there is only one possible velocity of propagation. As 
M. G. Kendall has remarked, ‘‘... the conclusion would be that the 
velocity of advance ultimately has the minimum value v,, whenever 
the initial distribution of mutant genes is strictly local in charac- 
ter.”? The rate at which the derivative of the initial distribution 
approaches zero as 2 approaches infinity is the critical factor de- 
termining the ultimate velocity of propagation of the wave of ad- 
vance. 

In the previous section we saw that our basic equation can be 
solved uniquely, provided that ¥,(z) exists and is L? integrable . 
Suppose that there were more than one velocity possible, then 
we would have a corresponding number of possible solutions to the 
dispersion-growth equation. If we assume that the specified con- 
ditions hold, then our uniqueness condition asserts that there is a 
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contradiction; hence there is only one possible velocity. This ve- 
locity is time-dependent during the transient state but should ap- 
proach a constant value as time passes (steady state). We are still 
faced with the question as to whether this wave front will be set 
up; again, with our uniqueness theorem holding, we can assert 
that there is some unique wave front in the steady state. This 
steady-state velocity should correspond to v,, but the author has 
not been able to prove this as a general statement. The stability 
of this wave front has already been shown by Skellam. If the 
initial-value function were such as not to be subsumed under our 
stated conditions, then we have no assurance that the velocity is 
v but may well differ from it (see Kendall, 1948, for a discussion 
of the analogous linear problem). Intuitively we can see what is 
happening by looking at (17). The A,’s are always positive, but, 
unless we have a bound on the order of the Fourier coefficients 
C,, we cannot assert that in the passage to the steady state the 
inhomogeneous term containing the initial values is insignificant. 

In this respect one has to reduce drastically the generalities by 
constructing concrete examples. In a paper in preparation, the 
Sturm-Liouville method is being applied to a number of interesting 
ecological situations. Unfortunately, a considerable body of ma- 
nipulation is necessary for quantitative results; however, for pur- 
poses of illustration we demonstrate the essentials of the method 
by only going to the second order of approximation. 

Illustrative example. We consider the case of a favorable habitat 
flanked by extremely unfavorable conditions at 0 and 7. Then we 
have to solve equation (1) subject to the boundary conditions 


w (0, ¢) = w(m, t) = 0 (¢2 0). (21) 


Let the initial-value function be a sine function: 
Ww vt)= Sin 2. 
° ( ) 


The diffusion and logistic growth coefficients are set equal to 
unity to simplify the arithmetic. 
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The normalized eigenfunctions ¢, (x) which satisfy the condi- 


tions (0) = d(7) = 0 are 


@, (2) = = sin iS a ( (ie wa SI 


The expansion of uw (za, ¢) is 


i oo 
ip («, t)~ — >3 U_(t) sin na, 


V7 


n=1 


as in (15). 
Now to expand f(z, ¢, w) = v’?- vu, 


f(2z,t,u) = = ye noe (t) sin NX, 


where 


F(é)= = - (u? — W) sin nada. 
et 


Substituting for the series (24) for w, we obtain 


P(t) =- 22 SOY) Ollamn) 0410) Ua (Os 


l=1 m=1 


here 


7 
Q(l,m,n) = ik sin Jz sin mx sin nad. 
i) 


(24) 


(25) 


(26) 


(27) 


(28) 
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Our infinite system of non-linear integral equations (17) now is 


—nt 


e Ar 
= e"* U_(s) ds - 
V7 Jo 


U,(t)= 5,7" + 


1 & = 
= j et) Y* Q(L,myn) U;(8) Up, (3) ds. (29) 


I,m=1 


We solve this infinite system by truncating to only the first three 


U_(¢). The truncated system is to be solved by successive ap- 
proximations. Letting 


U® (t) =0, n=1,2,3, (30) 
then the first approximations are, from (19), 


UD (t) = em', . 
31) 


US (t) = UM (t) = 0. 


The second approximation U(¢) is obtained from 


t qt) 
UM (= ayers | ernern Te a, 
0 7 


1 t 3 
5 i en) Y* Q(Lm,n) UP(s) UM(s) ds. (32) 


I,m=1 
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Carrying out the indicated operations, we have 
UP (t) = e“*(1 + 0740), 
Us” (zt) = 0, (33) 


US (t) = 


i (a3 = eo): 


Thus to the second order of approximation; from (24), we get 


ys? (x,t) = “: & ( + 2 sin x + = (en*h— ¢-**) sin 22].(98 


This approximation is good for only small values of ¢. Note that, 
when ¢ = 0, ¥™ (2,0) = J (2,0). Now W?(z,?) is a decreasing func- 
tion of the time (for fixed x); this can be interpreted as stating that 
dispersion is the dominating factor for small ¢ The reason is 
simple enough: the outflow due to diffusion just balances the lin- 
ear growth term, leaving the negative quadratic term W (z,{?). 
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An approximation method using a sine function is used to solve the 
second degree growth equation for the case in which an organism may 
simultaneously become dispersed throughout a uniform region. The re- 
sulting expression for a special case is compared with the expression ob- 
tained by R. Barakat (1959, Bull. Math, Biophysics, 21, 141-51), giving the 
first two terms, by an iterative procedure. The agreement is satisfactory. 


Recently, R. Barakat (1959) has considered the transient situa- 
tion for a growing, dispersing population. An example is treated in 
which there are absorbing barriers, the initial distribution being a 
sine function finite over the entire region. It is the purpose of this 
note to give an approximation solution for such situations using 4 
sine function instead of a linear function for the spatial distribu- 
tion of the population at any time. Non-linear approximations will 
also be used for other cases. In the notation used previously 
(Landahl, 1957), the differential equation for the number density (n) 
as a function of distance (x) and time (¢) may be written 


aman Bnd + D Vn, (1) 


One dimension. Absorbing barriers at +Z. Initial distribution, 
N=, COS 72/2Z. let us suppose that we may approximate the 
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solution for the population density n by the following expression: 
n(z,t) = m(t) cos 7 2/2Z, (2) 


a being a distance measured from the center of the region of width 
2Z, and m being a function of time. This expression satisfies the 
boundary conditions for n at the barriers, + Z, and is symmetric. It 
can be used as long as m(t) does not exceed n* = 4/B. The dif- 
ference, 


m(t + dt) cos 7 2/2Z — m(t) cos 7 x/2Z, 


integrated from —Z to Z represents the gain in numbers during the 
time dt. This gain must be equal to dt times the integral of an — 
Bn? from —Z to Z, the net birth rate minus death rate, minus 2D(dn/ 
dx)7dt, the loss due to the two absorbing barriers. Equating these 
two net gains and performing the operations indicated, we find, on 
re-arranging, 


; 7?D 7 
2 
n= (2) mE am (3) 


If the coefficient of m in expression (3) is zero or negative, then 
mis a decreasing function of time. Thus, m will decrease if 


Zea souhpelatier (4) 
“2 Va oo oe 


X being defined by (4). Solving (3) for this case and noting that 
m(0) = m9, the initial value of n at the origin, we find (Z < 1.11 X) 


—at 
ange@ cos72%/2Z z*D 
x,t) = ——-—______ "ge 


mm 1 —at 4Z? 
G+ 7rB nm -7 7B nye 


Xe (5) 


If Z < 1.11 X, a will be positive and the population will disap- 
pear. If the coefficient (—a) of m in (3) is positive, then m will in- 


POPULATION GROWTH 155 


crease according to (3). If this limit is less than or equal to n* = 
a/B, then the solution of (3), as given by (5), is a satisfactory 
solution. By setting m = 0 and m = n* in (3), we find, provided that 
Ny < n*, that (5) is a satisfactory solution if Z is not too large, 


D 
SAE RES ry Varn, = 240 x: 6 
(4 - w)a . 


Lee, 
From (5) it can be seen that n(0, ~), the steady state value of n 
at the center of the region, will be given by 


4 7 
m0, ~)-(=- = n* (7) 


in the range given by (6). In the case of the steady state it can be 
shown that, for Z less than about 1.8X, expression (5) gives a very 
close approximation to the exact solution (Skellam, 1951). 

If Z > 2.40 X, we use expression (5) for the solution until the 
time ¢* when n(z, t*) = n*. After ¢* let n(z, ¢) be given by 


nx, t)=n*, | «| <s 
eae ac (8) 


nz, t)=n* cos 72/(Z - 8), la|>s 


where s is a function of time determined by the condition that the 
integral of the net birth rate minus death rate over the region during 
dt minus that leaving at the absorbing barriers is equal to the 
change in total population given by subtracting the integral of 
n(x, t) from the corresponding integral in which s is replaced by 
3 +As. Performing these operations we obtain the following ex- 
pression for s: 


8, 


1 7 
(7 -2)8 = 5 0(Z- 8) (4-7) - 97 (9) 
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Since s = 0 for ¢ = ¢*, s is then given from 


2 2 
(Z-— 8)7-= Ee (22 Ea 


a a) eg 84> DG Tn 27 (10) 
a(4 - 7) a(4 — 7) 


Substituting Z — s from (10) into (8), we obtain the desired solution 
for ¢>¢*. Note that at t= 0, Z—~s =n VD/a(4 —7) =2.40X. 

We next compare the results of the above approximation method 
with that given by Barakat (1959) for a special case in which 
Ny /n* = 1/Vr, OD = 1, and. Z = at. Since a= 0 for this case, 
we obtain from (5) 


n(x, £) cos 2 


oS si 
n* Va +Ft a) 


Now in Barakat’s notation n(z, t)/n* is yw and cos z becomes sin 2, 
Since z is measured from one boundary instead of the center of the 
region as is done in (11). The time derivative of w or n/n* at the 
center of the region is —0.25 from (11) while the corresponding 
value from equation (34) of Barakat is —0.24. The coefficient of 
sin 3z in the latter expression is very small, so that the change in 
the form of the spatial distribution is negligible for small values of 
t. 

One dimensional infinite region. Sine approximation. Let the 
distribution be approximated by the following expression: 


1 1 
na, )/n* =~ cos n(x — s)/(r - 8)+55 (12) 


where r and s are functions of time to be determined. Initially r = 
™ and s=0. First we suppose that multiplication stops for a time 


A¢ and diffusion across the point z = 9" +s) moves r to r+Apr 
(cf. Landahl, 1957, Fig. 3). This requires that the integral from 


1 
5 + 8) to r+ Ar using (12) with r= 7+ Ar, minus the integral from 
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1 1 
5 (F + s) to r using (12), is equal to the diffusion flow across 5 + 


s) during A¢. Next, let diffusion stop but growth take place so 
that as a result s must change to s + As and the increase in the 
area under the curve is equal to the space integral of the growth 
rate times Az. These two equalities can be written as follows: 


(13) 


nd 
+8= i (r- 8). (14) 


Equations (13) and (14) are analogous to equations (21) and (20) of 
Landahl, 1957. [Note, however, that there is an error in equation 
(20) of loc. cit., a minus sign appearing instead of a plus sign on 
the left side of the equation.) Combining equations (13) and (14), 
we find the following expression which is analogous to equation 
(23) of loc. cit.: 


a PG ae ee eae 

= (r — 8) 5 tL (e—28)*], (15) 
rg _8a°D_ 1 
L SE ae (16) 


Solving (15) for (r—s) as a function of ¢ and substituting the re- 
sulting expression into (13) enables one to calculate r(¢). Thus 
with (r- s) and 7 known, we have s(¢). Hence the solution for 
nx, t) can be obtained from (12). Note that at ¢= 0, r—s=L’and 
since then 7 = 8, the population moves out with a velocity v = Hs) 
given by 


AEROS B 7 


irene kata eae (17) 
m-1L’ /8(7 - 1) 


) 
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a value rather close to that obtained using the linear approxima- 
tion. (Landahl, 1957, p.179). 

The above result suggests that one might use an exponential 
space distribution instead of a linear or sine function. Since n is 
non-zero everywhere, even without diffusion the population profile 
moves out as a wave with a velocity which approaches A, where 
\ is the space constant of the initial distribution. In order to use 
the approximation method with an exponential space distribution, 
the velocity obtained, which is the same as that using the linear 
approximation, must substantially exceed the velocity “A. Hence 
the space constant of the initial distribution, as well as that used 


: i 
for the approximate distribution must not be too large, i.e.,A < 3 bs 


If this condition can be satisfied, then the exponential approxima- 
tion can be used and the result is just that obtained by replacing 
the linear approximation at any time by a function which has cen- 


1 1 
tral symmetry about the point at hed + s) where n = 5% being made 
up of two parts. Each part is an exponential function of distance 
1 
having the same slope at n= 5" as the linear approximation, the 


part on the right approaching zero with distance while that on the 
left approaches unit. 

It is of interest to note the effect of the dividing point on the 
results obtained in the case of the linear approximation. In the 


. . 1 
previous cases we have divided the region at a de! Bee A= 


uJ 1 
5 ~ s). If we divide the region at r— 7" — s) we find that the 


limiting velocity in the linear case is VaD/3 instead of V2 «D/3. 


If instead of ~ op 2; ‘vidi 
instead oO g gin the above dividing point, we use 1/q, then 


the velocity approaches Vq%D/3. The transient expressions simi- 
larly will have modified time constants. 
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cussing the manuscript. 
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V. S. Ivlev [Experimental Ecology of Nutrition of Fishes, 1955, 
Moscow (in Russia)] has shown that the food uptake by fishes diving a 
fixed interval of time is an exponential function of the concentration of 
food. Ivlev’s equation is derived here, and it is shown that it can hold 
only for non-stationary conditions, such as prevailed in Ivlev’s experi- 
ments. For a stationary state, the rate of food uptake should tend asymp- 
totically to a limiting value as the concentration increases, but the vari- 
ation is not exponential. Different other aspects of the problem are in- 
vestigated, and definite new experimental procedures suggested. The 
implications of Ivlev’s findings on the effect of non-uniformity of food 
distribution upon the rate of food consumption are studied from a mathe- 
matical point of view. The conclusion is reached that whereas a fish 
does not, in the process of eating, move directly to an individual food 
particle which it perceives, it does move more or less directly to large 
aggregates of particles, if the latter are distributed nonuniformly. 


In his recent interesting book, V. S. Ivlev (1955) analyzes theo- 
retically and experimentally the relation between the rate of con- 
sumption of food by fishes and the concentration of the food. 

In the classical works of Alfred J. Lotka (1925) and of Vito 
Volterra (1931) on some aspects of interaction of different species, 
it is assumed that the total rate of change of available food is pro- 
portional to the product of the number of consuming individuals and 
the concentration of available food. Hence, per individual, the 
rate of consumption is assumed to be proportional to the concentra- 
tion of food. 

Ivlev (1955, pp. 19-20) correctly points out that such an assump- 
tion may hold at best as an approximation for small concentrations 


161 


162 N. RASHEVSKY 


of food. Taken literally, the assumption leads directly to a bio- 
logical absurdity, namely, the possibility of increasing the rate of 
food consumption of any individual beyond any limits. Actually, 
there is a limit determined by the physiological properties of the 
individual, such as rate of swallowing, delay of food in the stom- 
ach, etc. Denoting by r the rate of consumption of food, by p—its 
concentration, by 2—the maximum physiologically possible rate of 
consumption, and by €—a coefficient, Ivlev argues that the rate of 
increase of r with p, that is, dr/dp, must be proportional to the 
difference & —r. Hence the equation 


< = &(R -1). (1) 
P 


Integrating (1) we obtain: 
r=R(1-e-*"), (2) 


Ivlev has made numerous experiments to verify equation (2) and 
found it to agree amazingly well with his experimental data (1955, 
p. 21, Figure 1). From the experimental 7,p curves Ivlev calculates 
the coefficient € for different cases. The physical meaning of & 
remains, however, entirely undefined because of the purely formal 
nature of the mathematical approach. 

One of the purposes of this note is to analyze the possible bio- 
physical meaning of equation (2). 

The first thing which suggests itself is to connect the rate of 
eating with the speed of locomotion of the fish. A hungry animal 
is likely to move towards a food particle which it sees much more 
rapidly than a well-fed animal. The distance between the food 
particles decreases with increasing concentration; therefore, for a 
constant speed of locomotion, the interval between the swallowing 
of food particles will decrease with increasing concentration, and 
the rate of eating will at first increase with p. However, as the 
animal eats faster and faster, it rapidly becomes satiated and its 
speed of locomotion decreases, thus increasing the intervals be- 
tween the swallowing of food particles. Expressions which show 


an asymptotic approach of r to R can be derived from a further 
elaboration of this picture. 
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Unfortunately the picture does not correspond to actual observa- 
tion (Ivlev, 1956). A hungry fish, as a rule, does not move faster 
than a well-fed one. Moreover, a fish never moves from one food 
particle to another in a straight line, but stops, makes turns, 
passes some food particles, etc. Professor Ivlev (1956) describes 
the behavior of fishes as erratic and capricious. It is this very 
feature that may give us a clue to the solution of the problem. A 
probabilistic approach is definitely suggested. 

Let us assume that the fish moves in a random manner, with a 
constant average velocity v. Let those movements have nothing to 
do with the distribution of food particles. We may possibly view 
the situation in the following manner. A fish moves in a random 
fashion amongst uniformly distributed particles. Occasionally a 
food particle comes within the animal’s range of perception, and it 
may eat the particle. The swallowing and the digestion of the 
particle take a finite time, 7; therefore, we may assume that if, 
after swallowing a particle, the fish finds another one after a time 
less than r, it will not eat it. If the concentration of the particles 
is so small, that the interval of time between two particles coming 
within the range of perception of the fish is appreciably larger than 
r, then every particle perceived will be eaten, and the rate of con- 
sumption of particles will be proportional to the concentration of 
the particles. If, however, the concentration of the particles is 
very large, then the fish will not swallow every particle it per- 
ceives. The best it can do is to swallow a particle every; sec- 
onds. Hence, no matter how much the concentration of the parti- 
cles increases, the fish can swallow at most 1/7 particles per 
second. 

We shall find subsequently the type of mathematical relations 
between concentrations of particles and the rate of their consump- 
tion, to which the above picture leads. First, however, we must 
discuss another important point. The assumption of proportion- 
ality made by A. Lotka and V. Volterra contains tacitly another 
assumption: namely, that the whole process of food consumption is 
stationary, and that the conditions of the consuming organism do 
not vary during the process. If those conditions should vary, this 
would be reflected in a variation of the coefficient of proportional- 
ity between rate of consumption and concentration of food, and 
thus would in effect lead to deviation from exact proportionality. 
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The same tacit assumption is also quite compatible with Ivlev’s 
argument which leads to equation (2). The assumption of a station- 
ary regime also would underlay the derivation of any equations 
based on the picture which we just outlined. 

If, however, we read carefully the description of Ivlev’s experi- 
ments, we find that they did not occur under the conditions of a 
stationary regime, and that the latter was apparently not even in- 
tended. Small fishes were first kept without food for 18-20 hours 
(Ivlev, 1955, p. 21). After that 5 fishes were placed in a shallow 
tray in which food particles were distributed uniformly. The num- 
ber of particles swallowed by each fish during a period of 1% to 2 
hours was directly counted (loc. cit., p. 20). Also, after the end of 
the experiments, that is after 1% to 2 hours, the content of the gut 
of each fish was determined by autopsy. 

The experiments were conducted for different concentrations of 
food particles. The concentrations varied within a ratio from 1 to 
10. The amount of food consumed during the time of the experi- 
ment was plotted against the concentrations, and an excellent 
agreement with equation (2) was found. 

Three different types of fishes were used, their average weights 
being correspondingly 1349 mg, 1826 mg, and 673 mg. The maxi- 
mum (asymptotic) amounts of food consumed by the three types of 
fishes were correspondingly 292 mg, 198 mg, and 90.0 mg. 

Comparing the weight of the food inside the fish at the end of 
the experiment with the weight of the corresponding fish itself, we 
may safely conclude that towards the end of the experiment the gut 
of the fish must have been almost, if not totally, filled with food. 
On the other hand, because of the preceding period of starvation, 
at the beginning of the experiment the guts were empty. The feel- 
ing of hunger must certainly have been much greater at the begin- 
ning of the experiment than at the end. Therefore, it is only natu- 
ral to assume that for the same concentration of food, the fish 
might have eaten much faster at the beginning than at the end. 

Let us see how the situation is described mathematically. 

Let a food particle become noticeable to the fish, either through 
sight or through olfactory sense, when it is within a distance o 
from the fish. If p, the concentration of food, is expressed as the 
number of food particles per unit volume, then the number of food 
particles which will come to the attention of the fish per second is 
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given by 
N= 71VO0'D. (3) 
We thus have 
n= AD, (4) 
with 
Q=1o0'v. (5) 


In seeing or sensing a food particle, the fish may or may not eat it, 
depending on the intensity of its hunger. Somewhat in analogy 
with the accepted (though by no means completely proved) notion 
about the nature of hunger in man, we may quite plausibly assume, 
that the emptier the gut of the fish, the more hungry it will feel, 
and the greater therefore the probability of its swallowing a food 
particle which it senses. If M denotes the maximal possible con- 
tent of the gut in food, if m; denotes the actual amount of food in 
the gut, and if a@ is a factor of proportionality, then the probability 
of the fish eating a particle ‘‘at sight’? is equal to a(M — m,). When 
m,=0, that is when the gut is completely empty, we may reason- 
ably assume that the fish will eat any food particle which it 
senses. Hence the probability must be 1 for m;=0, or a=1/M. 
Thus the probability P of the fish to eat the particle which it 
senses is: 


M—™m™,; 
M 


P= (6) 


The fish senses n particles per second, where n is given by (4); 
therefore, the total number of particles eaten per second is equal 
ton=nP =aPp. Let m, denote the mass of a particle. Then the 
total amount 7 of food ingested per second is equal, according to 


(5) and (6), to: 
apm , 


[= TaTEE (M - ™M ;)- (7) 
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The content of the gut is emptied continuously by absorption of 
digested food, and discontinuously through defecation. From the 
circumstance that Ivlev used as a control the direct determination 
of food in the gut after a 2-hour long experiment, it follows that 
during this time no appreciable emptying of the gut did occur. In 
other words, during the experiment the gut was merely being filled 
through ingestion. Therefore, with the same degree of approxima- 
tion as we can neglect the emptying of the gut, we may put: 


dm, dpm, 
dt M 


(M-m,). (8) 


Integrated with the initial conditions m, = 0 at ¢ = 0, this gives: 


ampP 
es Bae 
m,=M(1-e ). (9) 


For a constant p, the quantity m,; is a function of ¢, which ap- 
proaches asymptotically M. But for a constant t, m; is a function 
of p, which approaches the same asymptote. Ivlev kept the dura- 
tion of his experiments constant, about 144-2 hours. Let that 
duration be denoted by ¢,. He furthermore identified the amount Mm; 
consumed during the constant time ¢, with the rate of consumption, 
which was tacitly assumed constant during the time ¢. Hence 
Ivlev’s ris actually our m;. Putting 


amt 
m,=t; M=R; Pp 


=, (10) 


we obtain from (9) Ivlev’s equation (2). 

According to the above, if Ivlev would have varied ¢, within a 
wide range, say from 1 to 5 hours, and studied the variation of his 
curves with ¢, as parameter, he would have found different values 
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of € for different values of ¢,, but his exponential relation would 
still hold. Actually, what Ivlev measured is not the instantaneous 
rate of consumption /, as given by (7), but the average rate during 
the period ¢,, during which the instantaneous rate varies as 


_ &mpp : 


Ampe MM . (11) 


which is obtained by differentiating the right side of (9) with re- 
spect to ?. 

If our interpretation of Ivlev’s result is correct, then we should 
be able to estimate at least the order of magnitude of the constant 
€ in Ivlev’s experiments. Since no adequate data which we need 
for this purpose are given by Ivlev, we must restrict ourselves to 
indirect estimates. 

Strictly speaking, the above theory applies only to Ivlev’s third 
series of experiments, in which Daphnia pulex was used as a prey. 
Though Ivlev used shallow trays, yet the Daphnia would be dis- 
tributed throughout the whole volume of water, and the ‘‘collision’’ 
expressions (4) and (5), which are patterned after the kinetic the- 
ory of gases, hold. In the first series of experiments denatured 
fish eggs were used as a prey. In the second series live larvae of 
Chironomidae were used. In both cases the prey would be located 
at the bottom of the tray. The fish has to ‘‘dive’’ in order to pick 
it up, and relations (4), and especially (5), would not hold. We 
shall see subsequently how we may try to circumvent this diffi- 
culty at least in a very rough manner. Now let us consider the 
third series of experiments. 

It seems reasonable to assume that for fishes of the size used 
(total average weight 673 mg), the value of o is of the order of 
fraction of a cm. If for the velocity v we take the values from a 
previous paper of Ivlev (1944), we may put v ~ 1 cm sec™*. This 
gives & ~ 107* cm® in absolute units. But Ivlev gives his con- 
centration not in numbers of particles per cm’, but in milli- 
grams per cm’ of the tray. Nothing is said about the size of 
the tray. However, Professor Ivlev informs the author that the 
trays were about 50 x 60 cm, the water being about 4 cm deep. The 
weight of a Daphnia may vary between a few milligrams to, say, 
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10 mg. Let us take m,~1 mg. Then 1 mg cm~* actually meant 
0.25 mg cm7* of water volume, and in terms of number of Daphnia, 
it meant about 0.2 food particles per cm*’. Hence, to give the cor- 
rect number of collisions our & must be multiplied by about a factor 
of 0.2 making « ~ 0.02. As M we may safely take the maximum 
amount of food consumed by the fish. Thus we have: « ~ 2 x 107%; 
m, ~ 10° gm; M~ 0.1 gm; t, ~ 7,000 sec. Equation (10) now 
gives ¢~ 1.4, or recomputed to decimal exponential, as Ivlev uses, 
€~ 0.6. This is 3 times larger than Ivlev’s value €= 0.1917. If 
we increase &, by assuming o and v to be larger, the discrepancy 
becomes worse. The quantities o and v are the most dubious ones 
in the case.’ On the other hand, if as Professor Ivlev informs the 
author, a fish has sometimes to actually ‘“‘bump into the food,”® and 
if this holds for this particular series of experiments, then o would 
be reduced to about 0.1 cm, o? is reduced to about 107?, and the 
corrected value of « (reduced to concentrations expressed in milli- 
grams) becomes about 4 x 107*. This will give € in Ivlev’s units, 
as 0.12, which is not bad at all. 

The result is not sensitive to Ty s because the corrected value of 
& is obtained from the one expressed in absolute units by dividing 
the latter by 4 m,, 4 cm being the assumed depth of the tray. 

It may also be questioned whether the velocity of movement of 
the fishes in experiments of the third series is even approximately 
the same as the velocity of movement of the fishes used by Ivlev 
previously (1944). It must also be noted that, when comparable 
food was used (Cyclops), the unit concentration corresponded to 
40 specimens per liter of water (Ivlev, 1944, p. 139), while in the 
present case it appears to have been as much as one specimen per 
4 cm® or 250 per liter. 

The above discussion illustrates the urgent need of determining 
and keeping careful track of a number of biophysical parameters, 
which were completely neglected by Ivlev because of his purely 
formal approach. The above mentioned possibility of even a rough 
agreement between calculated and observed values of é indicates 
in our opinion that the. mathematical biophysical approach is not as 
hopelessly more difficult than the formal one as it may appear at 
first sight. ' 

In dealing with fishes which pick up their food from the bottom, 
we cannot use expressions (4) and (5), but must use a different 
approach. If we stick to our assumption that the movement of the 
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fish is not directed but random, we should consider the movement 
of the fish as a Brownian movement of a particle. It should be pos- 
sible to determine experimentally the average length of run between 
two turns and from this to compute the average displacement of the 
fish in a given direction. Then we may compute the average number 
of downward movements, in which the fish would strike the bottom 
with its nose. Knowing the size and number of food particles per 
cm’, we can then calculate the average number of collisions per 
second. It will be again of the form (4), that is: 


A Nil Bp, (12) 


and equation (2) will still hold for not too small concentrations, 
but 8 will be quite different from « in (4) or (5). 

In a shallow tray, however, the fish is likely to have a greater 
tendency to move horizontally than vertically in a random motion. 
We thus have a combined effect of randomness and directedness. 
The interesting discussion by N. J. Kobozev (1948) on vectorial- 
Brownian processes may be pertinent to such a study (cf. also 
Patlak, 1953). 

A very rough estimate of the value of € in the first two series of 
Ivlev’s experiments, in which the fish had to dive to the bottom for 
food, may be made as follows: 

Assuming that the fish moves with a velocity v with equal proba- 
bility in all directions, we may roughly consider that it moves with 
a velocity v/3 in the vertical direction. With v ~ 1-3 cm/sec, and 
the depth of the tray ~ 4 cm, the fish will move up and down, mak- 
ing a total of 8 cm, during about 10 seconds. Thus it will hit the 
ground once in every 10 seconds. In the first experiment, fish eggs 
were spread with a lowest concentration of 1 mg cm’ of tray. Each 
egg weighs according to Ivlev about 10~* mg. Hence, at the lowest 
concentration used, there are about 10 eggs per cm’. The cross 
section area of an egg is about 0.2 mm? or 2 x 107° cm’. If the fish 
has to bump with its nose not more than 1mm from the egg in order 
to find it, then the effective ‘‘target area’’ for each egg is about 
4 x 107? cm?. There are 10 eggs per cm’; therefore, the probability 
of the fish hitting an egg at each dive is ~ 4x 10-*. There is one 
dive per 10 seconds and therefore altogether, 4 x 10-? collisions 
per second. By definition B (equation 12) is essentially the num- 
ber of collisions per second for a unit concentration. But Ivlev 
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used 1 mg cm7? as unit concentration. Hence B ~ 4 x 10~*. With 
M~ 0.3 gm., m, ~ 10-* gm, ¢, = 7x 10° sec., we find from (10) in 
which 6 is substituted for 4, €~ 0.1, or in decimal exponents, 
& ~ 0.2, which compares favorably with Ivlev’s € = 0.16. 

A similar consideration for the second series of experiments 
leads to similar results. All we can expect is agreement of orders 
of magnitude. 

So much for the non-stationary regime. What happens, if the ex- 
periment is carried out for an indefinite time, when a stationary 
regime is reached? In a stationary regime the maximum rate of con- 
sumption of food is determined by the rate of emptying of the gas- 
trointestinal tract either by absorption or by defecation. Neglect- 
ing the latter, we may put the rate of the former as equal to km,, 
where & is a constant, that can be determined experimentally from 
purely physiological experiments. In that case, we cannot neglect 
km,, and instead of (9) we must write: 


dm, am, Pp 


ee Oe (M —m,) -km,, (13) 
or 
dm; Am , p 
“Betis Piet iu +k] m;. (14) 
Integrated this gives: 
am P -( to k)t 
iS aaa = ). (15) 
Raed Li ey 
m 


Since during 2 hours, or about 10* seconds, there was no ap- 
preciable emptying of the gut, we may conclude that k << 107+ 
sec and is not greater than 10-* sec~'. As we have seen, at the 
very lowest concentrations n~ 10-' particles per cm, or less. 
Even with the smallest value of & chosen, namely « ~ 107? in ab- 
am,” 


solute units, we find that , even for the smallest concentra- 
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tions used, is about 10 & For the more plausible higher values of 
a, the inequality 


>>k (16) 


is even larger. Thus equation (15) reduces to (8). If Ivlev had 
gone to much lower concentrations, for which (16) would not hold, 
he should have found appreciable deviations from the exponential 
relation (2). 

However, no matter how small &, the gut is never filled com- 
pletely. The filling approaches asymptotically the value 


a, (17) 


which may be only one food particle less than the full value M, but 
it zs slightly less. Now we must remember that in equation (14) 
dm ,/dt does not any more represent the food intake. It is the dif- 
ference between the intake /, given by (7), or by the first term of 
the right hand side of (13), and the loss km,. The rate of intake is 
still given by (7). 

After an infinite lapse of time ¢ a stationary state is reached, in 
which dm ,/dt = 0, or 


am, p 
M 


(M -—m,) = km; (18) 


But in this state m; is given by (17). Hence the rate of consump- 
tion (7) is given in the stationary state by: 


am, DEM 


fn ee Lee (19) 
rot am, p + kM 
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For very small values of p, /,,,, = %m,p, it is proportional to 
the concentration p. As p — ~, /.,,, —> kM, which is the rate of 
absorption from the gut as should be physically the case (when 
defecation is neglected). 

Thus for a stationary regime the Lotka-Volterra assumption holds 
only for small concentrations of food. With increasing concentra- 
tion the rate of absorption tends to an asymptotic value which is 
determined by the physiological constants of the organism. The 
range of concentrations in which this variation of rate of consump- 
tion is noticeable, is, however, for stationary states, much below 
the concentrations used by Ivlev. 

An equation of the form (19) can be derived for the stationary 
state also by a different consideration. We shall follow an argu- 
ment used by John Calhoun (1957) in an entirely different connec- 
tion. 

Consider a sufficiently long period of time ¢. Let n’ be the num- 
ber of particles eaten by the fish per unit time. Then, during the 
time ¢ the fish will eat n’¢ particles. If after eating a particle the 
fish does not eat for a period 7, then during the time ¢ it will not 
eat any particles for a total period of n’tr, regardless of whether 
during this period it encounters any particles or not. It will there- 
fore eat all particles it encounters during the time ¢-n’tr=¢(1- 
n’r). Per unit time it encounters n particles; therefore, during the 
time ¢(1-—7’r) it will encounter altogether nt(1-—-n’r) particles 
and, as just said, it will eat them all. Hence during the time ¢ it 
will eat altogether né(1—-—n’r) particles. But, as we remarked 
above, during the time ¢ it eats n’t particles. Hence 


nt (1 —n’r) = n't, (20) 
or 


n(l-n’r)=n’'. (21) 


Solved for n’, equation (21) gives 


n 


, 


n= 


1+ nz co 
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Introducing (4) into (22) we find: 
If m, is the mass of the food particle, then the amount r eaten is 
given by: 


om , P 


1+orp u») 


Equation (23) shows that for very small values of the concentra- 
tion p of food, the rate of consumption r increases as am, p, that is 
proportional to the concentration. As p increases, however, 7 tends 
asymptotically to the value: 


ep 
ro= —, 
oe. (24) 
If we put 
m 
Pop: ae 
ona am, = €h; (25) 
then equation (23) becomes: 
R 
1+ €p 


Equation (26) is basically of the same form as (19). 

For small values of p, r increases as &p, just as it does accord- 
ing to equation (2), and for very large values of p, the rate of con- 
sumption tends to R, as it does according to equation (2). How- 
ever, while the general shape of the curves represented by equa- 
tions (2) and (26) or (19) are similar, equations (19) and (26) do not 
represent quantitatively Ivlev’s experimental data. The curve of 
equations (19) and (26) runs much lower than the curve of equation 
(2), for the same initial slope and the same asymptote. According 
to (2), the asymptotic value F is practically reached for much 
smaller values of p, than according to (26). 
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All of this is of course not surprising, because equation (2) or 
its equivalent (9) represent the state of affairs for a nonstationary 
regime, while equations (19) or (26) represent the situation for a 
completely stationary regime. To the extent that Ivlev’s experi- 
ments were carried out under non-stationary conditions, there is 
thus far no evidence whatsoever that (19) or (26) contradict any 
actual facts. Only further experimentation can decide this. 

It may, however, be of interest to see what consequences will 
have to he drawn if further experiments reveal that an equation of 
the form of (2) holds also for stationary regimes. 

In deriving equation (26) we assumed that the quantities « and 
7 are constant. The constancy of & implies either the constancy of 
o and of v, according to (5), or else it implies an inverse propor- 
tionality between o” and v. The latter assumption is obviously ab- 
surd, and we may dismiss it. The assumption of a constant v is 
borne out by earlier experiments of Ivlev (1944). 

In an earlier paper, Ivlev (1944) studied the relation of the speed 
of locomotion of the fish to the concentration of food. He found 
that this speed is essentially constant, except for very high con- 
centrations. In those experiments, Ivlev used concentrations which 
varied from 40 to 1000 food particles per liter (see above, p. 168). 
In his later experiments, made for the purpose of verifying (2), Ivlev 
used concentrations from 250 to 2500 food particles per liter. 
Therefore, the assumption of a constant v in our derivation of (2) 
is not entirely justified, because v begins to decrease appreciably 
already for concentrations above 600 particles per liter whereas 
the asymptotic level R is reached for concentrations around 4000 
particles per liter. However, in the two cases different species of 
fishes and different food particles were used. Therefore, roughly, 
the assumption of a constant velocity may be justified. 

Hence there remains only the possibility of a variation of o or of 
r, or of both. A variation of o is not very likely. We shall return 
to this later, and eliminate the possibility entirely. Now we shall 
discuss a possible variation of r. 

It would be rather natural to expect that 7 will depend on the rate 
r of consumption of food. For a given average rate of emptying of 
the stomach into the intestinal tract, the stomach will be the fuller, 
the greater r, The more full the stomach, the more likely that a 
fish will wait for a longer time between the eating of two particles 
of food. Hence we might expect that 7 will increase with r. This 
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is indirectly borne out by Ivlev’s earlier observations (1944), that 
a fish takes a longer time to swallow a food particle when the con- 
centration of food particles is large. We shall determine the shape 
of r(r) which would lead tothe equation (2) instead of equation (26) 
which is obtained for a constant r. 

From equation (25) we notice that for a constant 7 the maximum 
obtainable rate R of consumption is simply 1/7. Intuitively this is 
quite clear: no matter how large the concentration p of food, the 
fish cannot eat more than 1/r particles per unit time. 

The unit of time to which the measured rate of consumption is 
referred must of course be chosen sufficiently large, in order that 
N = 1/r would be large. This is also necessary from the purely 
experimental point of view, regardless of any theory. The unit of 
time must be large enough for the fish to consume, even at the low 
rates, a large enough number of food particles. 

Now, however, that 7 is itself a function of 7, we cannot use the 
expression (25) R = m,/T. However, since FR € in (2) and Om. in 
(23) are both equal to the initial slope of the curve, therefore we 
can still put am, = ER, 

Since 7 is not now a constant, but a function r(r) of r, weno 
longer obtain equation (26) from (23). Equation (23) solved for p 
gives: 


pos: 
p= aN. (27) 


If equation (23) for a proper choice of r(r) should lead to (2), 
then r(r) must be chosen so as to make the right sides of (27)and 
(28) equal. Hence r(r) is determined by the equation: 


if R 
db. io. 2 ie = log ——. (29) 
Kobrin .€ R-f 
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Remembering that am, = ER, according to (25), and solving for 


r(r), we find: 


By ls ad SOU A Sa oe (30) 
: 


= 
R tog (1 - =) 


For very small values of r/R we may expand the logarithm in the 
denominator, breaking up after the quadratic term. We then find 
after simple rearrangements: 


7(r) = ——_——_.. (31) 


As r tends to zero, this expression tends to 


ate 


0) = 32 
(0) = <2 (32) 
On the other hand when r = RF, equation (30) gives: 
(R) = = (33) 
Tt — R . 


Thus r(r) increases from m,/2R to m,/R as r increases from 
zero to its maximal possible value R. In other words PF is the in- 
verse of the maximum value of r. 

Such a variation of r with r does not appear to be implausible. 
However, the circumstance that r(r) changes exactly by a factor of 
2 between r = 0 tor = R, appears to be rather strange. In principle, 
equation (30) could he verified experimentally in the following 
manner. 

A fish, confined to a small space, is presented food particles at 
a given rate r per second. The rate is chosen so as to be below 
the maximal rate of consumption of the fish. In other words, the 
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rate should be small enough to make sure that the fish eats every 
particle presented. After a stationary regime is established, the 
fish is presented occasionally with a food particle after a time 
t* < 1/r, following a swallowed particle. In this manner it may be 
possible to establish that for a given 7, the fish will not eat a 
particle if presented after a sufficiently small time ¢*. Fora given 
r, the smallest time ¢* for which the fish will still eat the particle, 
is the corresponding value of r(r). The experiment may not be 
practical. It seems more likely that a somewhat less direct but 
easier verification of (30) could be obtained by taking motion pic- 
tures of the fish during the process of feeding and by analyzing 
then the detailed mechanism of approach, swallowing, or rejection 
of a food particle. 

We mentioned earlier the possibility of a variation of o. Though 
unlikely, it may be thought that with a very large number of parti- 
cles, the fish fails to see all of them, being, so to say, distracted 
by too many of them. The phenomenon would have a neurobiophysi- 
cal basis in N. Rashevsky’s theory of mutual inhibition of too 
many stimuli (Rashevsky, 1959). If such a thing should happen, 
we would expect a decrease of o with n. This, because of (5), 
would mean a decrease of & with p in equation (33), or a decrease 
of € with p in equation (26). If, however, we put € = €(p) into (26) 
and determine by the similar argument as was used for 7, the shape 
of the function &(p) which would make (26) to acquire the form of 
(2), we find that € actually must increase with p. This is an ob- 
viously absurd situation, which eliminates the possibility of a 
variation of o. 

All our discussions hitherto were based on the assumption that 
the fish picks up its food during a process of random movements. 
Though this is borne out by observations, and according to Ivlev, a 
fish never moves directly towards a food particle, it may be of in- 
terest to see what type of relations are obtained if a directed move- 
ment is assumed. We shall consider just two examples, 

Let the average distance between two food particles be /. The 
fish swallows a particle and moves to another with a velocity v. 
We may also assume for generality that the hungrier the animal, the 
greater the velocity with which it moves toward the particle. Then 
the velocity v is given by 


v= v,(M-™m,), (34) 
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where v, is the velocity of a completely starved fish. 
Per unit time the fish will eat v/l particles. Hence for a non- 
stationary case, as exemplified in Ivlev’s experiments: 


ere eg he (35) 


If p is the number of particles per cm* then we have approxi- 
mately: 


Eee. (36) 


Introducing (36) into (35) and integrating with initial condition 
m,=0 at ¢=0, we find: 


BPS Wer a (37) 


If the experiment is carried out for a fixed time ¢ = t,, and if we 
put 


m,=?t;M=R; m,%,t, = €, (38) 


we find 


r= R(1-27* ’?), (39) 


This is different from (2) and will not agree with Ivlev’s data. 

For a stationary regime, we may argue as follows: in this state 
the gut of the animal is almost full, and the fish will run after 
another food particle only after a time r following the swallowing 
of food. Because of its satiation, the velocity » of the animal will 
also he practically constant. Hence the average time which it will 
take between the consumption of two food particles is the sum of 
r and of the time //v which it takes to move to the next particle. 
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Hence per unit time the number of particles swallowed is 


or, according to (36) 


vp 
Te ey) 


The average amount r of food consumed per unit time is 


m v%/p 


r= ———__,—.. 41 
Ll+rv%/p oy 


For small values of p, r varies as %/p. The tangent of the z, 
Mm: 
p-curve is vertical at the origin. As p > ~,r > —. 
fi 
In his book Ivlev also studies the dependence of the rate of eat- 
ing on the nonuniformity of the distribution of the food. Food parti- 
cles were arranged in different patterns and as a measure of non- 


uniformity the quantity 


c= // = (42) 


was used. Here £ denotes the absolute value of the deviation 
from the average concentration in a given region and n—the num- 
ber of regions. By a somewhat formal mathematical argument simi- 
lar to the one which led him to (2), Ivlev arrives at the following 
expression: 


p= R [1 — e7& P+ 5], (43) 
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which, he also finds, represents well the experimental data. For 
€ = 0, (48) reduces to (2). 

At first we may be tempted to derive equation (43) from (2) by 
considering the average effect of a nonuniform distribution. It is, 
however, readily seen that in such a case r should decrease with 
increasing ¢, whereas it actually increases. 

Let 7 be expressed quite generally as a function of p by the 
equation: 


r= f(p). (44) 


Let p be a function of the coordinates z and y. (Two-dimensional 
situations were used in Ivlev’s experiments.) 


p=u(a, y). (45) 
Kiquations (44) and (45) give: 
r= flu(2, y)). (46) 
Let p, be the average concentration, and put 
U(z, Y) — Po = u,(#, y)- (47) 


We then have: 


Jf (x, y) dx dy =0, (48) 


where the integration is extended over the whole experimental 
region. 
Consider now the average value 7 of r: 


1 
F= fre, y)] dx dy, (49) 


where S is the total area of the region over which the integration is 
extended. Introducing (47) into (49) we find: 


1 
f= {rt + U, (x, y)] de dy. (50) 
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Expanding f into a Taylor’s series, we have: 


et 1 (¢(df 
rm =| [tive de dy + Se u, (2, y) de dy + 


1 df 
aie u2 (x, y) da dy +... (51) 


Here the subscript p, denotes that the values of df/dp and d?f/dp* 
are taken for p=p,. The first integral of (19) equals f(p,). We 
thus have 


* 1 /d 
F=f (Po) + 3 () Jf (w, y) dx dy + 


Re Bie 3 
5 a) flew Yala dy tres (52) 


Because of (48), the first integral in (52) vanishes. The second 
is essentially equal to €?, as defined by (42). Hence, 


_ 1 fd ; 
F=f (Po) + 5 i €? +... + higher terms. (53) 


But if f is given by either (2), or (19), or (26), or by any other 
similar expression, then d?f/dp? <0, and therefore 7 decreases 
with increasing ¢, contrary to observation. From (53) it is also 
seen that in general 7 depends also on the higher moments of the 
distribution of particles. For symmetric distribution only odd mo- 
ments remain. 

The conclusion seems inevitable that the increase of 7 with fis 
due to the fact that the fish selects regions of highest food concen- 
tration, as Ivlev suggests (1955, p. 24). The maximum concentra- 
tion will be in general proportional, at least for not too unusual 
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distributions, to €. Thus 
Pmax = Pork. (54) 


In expression (2) we must substitute p_., for p. We then obtain 
an equation of the form of (43), with &€ = x. 

Thus though the fish does not move directly towards a food parti- 
cle when the latter are distributed uniformly, it does move towards 
a larger aggregate of such particles, if this aggregate is character- 
ized by a higher concentration. If the principal stimulus to the 
fish is a visual one, this may be readily understood: the sight of a 
large mass of particles produces a stronger stimulation than the 
sight of a single particle. If the stimulus is olfactory, the phe- 
nomenon can also be understood. A large number of quasi-uniformly 
scattered particles produces only small local olfactory gradients to 
which the animal may not respond. A large cluster of particles 
produces a gradient of a different order of magnitude. Such a 
gradient may well evoke tropistic responses. 

The author wishes to express his special thanks to Professor 
V. S. Ivlev for continued discussion and many valuable suggestions 
in the course of this investigation. The author is also indebted to 
Professor Herbert D. Landahl for a critical reading and discussion 
of the paper. 
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Some relational aspects of the property of self-reproduction of biologi- 
cal systems are studied. If in addition to the requirement of the property 
of self-reproduction we add also the requirement of adaptability of the 
organism to changing environment, this imposes certain conditions on 
the topology of the graphs which represent such systems. A further 
study of the relational properties of such systems seems to offer the 
possibility of deriving the principle of biological mapping from the re- 
quirement of self-reproduction and adaptability. 

An examination of the problem of the original formation of such self- 
reproducing systems in connection with the established fact of impossi- 
bility of spontaneous generation leads to the conclusion that an organ- 
ism must inhibit such processes which, in the absence of organisms, 
would lead to spontaneous generation. 


One of the most important characteristics of life is the property 
of every living organism to multiply by self-reproduction. Another 
equally important characteristic is that self-reproduction is the 
only method of multiplication thus far known. Multiplication of 
plants from cut-off parts is still fundamentally due to the self- 
reproduction of the cells which constitute these parts. Omnis 
vivum e vivo appears to be as true now as ever. 

There are two schools of thought in regard to the property of 
self-reproduction. One considers self-reproduction as the property 
of the whole complex physicochemical system which constitutes 
even the simplest organism. Self-reproduction, according to this 
point of view, is essentially due to the interactions of different 
parts of the complex biological system rather than a characteristic 
of one particular part of it. The other point of view, recently de- 
veloped mainly as the result of study of viruses and of some ge- 
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netic phenomena, ascribes the property of self-reproduction to a 
particularly structured ‘‘living molecule’’ which has the property 
of replicating itself. 

Our earlier studies on various possible causes of division of 
metabolising systems (Rashevsky, 1959) based on diffusion drag 
forces or any other similar concepts, are essentially representa- 
tives of the ‘‘systemic’’ approach to the problem of multiplication. 
Thus far, however, those attempts met only with a very limited 
success (Rashevsky, 1959a), 

In his well-known book A. I. Oparin (1957) adduces numerous 
convincing evidences for the systemic point of view, though this 
point of view remains still far from being established. Neither the 
proponents of the systemic point of view, nor those of the “tmolec- 
ular’’ one, have actually suggested any specific mechanism of 
self-reproduction, nor studied the possible general properties of 
such mechanisms. A self-replication of a molecule is supposed to 
be based on the ‘‘template principle.’? But why only extremely 
large molecules can act as templates and why we do not find self- 
replication of that type in simple molecules are questions that re- 
main unanswered. N, Rashevsky’s (1959c) recent discussion of 
some possible mechanisms of such a molecular self-reproduction 
strongly suggests that a self-reproducing molecule would have 
that property only in the presence of some definite structures, the 
formation of which must itself depend on the presence of the self- 
reproducing molecule. Thus we practically are led to the sys- 
temic approach. Even if the self-reproduction of a molecule takes 
place without that molecule forming a part of a larger system, the 
very complexity of such a self-reproducing molecule and the very 
large number of atoms of which it must be composed would sug- 
gest that the molecule itself must be considered as a complex 
system and that the process of self-reproduction is the result of 
an interaction of the different parts of this “*single-molecule”’ 
system, 

In the present note we shall discuss some relational aspects 
of self-reproduction from the systemic point of view. This, how- 
ever, Should by no means be construed as a rejection on our part 
of the molecular point of view. 

In a previous paper (Rashevsky, 1954) we introduced a repre- 
sentation of an organism by an oriented or directed graph. The 
vertices of such a graph represent the different biological proper- 
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ties of the organism, while the oriented edges of the graph repre- 
sent the relation of immediate precedence or immediate succes- 
sion between a pair of biological properties (Rashevsky, 1958). 
The biological properties are ‘‘carried’’ by some specific parts of 
the organism, those parts being either individual cell-types or 
whole organs (Rashevsky, 1954, 1956b, c). 

A somewhat different representation of an organism by a graph, 
a representation which in some respects proves to be more advan- 
tageous, has been recently introduced by Robert Rosen (1958a). 
Rosen considers the vertices of the graph as representing directly 
the different parts, or components, of the organism, while the con- 
necting lines represent different biological inputs and outputs into 
and from those parts. In many respects the two views are inter- 
changeable and either one can be used. It is more convenient 
here to use Rosen’s point of view, and we shall do so even when 
referring to conclusions which we made originally on the basis of 
our concept of the organismic graph. 

As we pointed out (Rashevsky, 1954, pp. 326-27), the exist- 
ence of directed (Rashevsky, 1955c) cycles in the graph of an or- 
ganism reflects the impossibility of an organism to form spontane- 
ously without another similar organism already being present. 
From a purely formal mathematical point of view a graph with a 
directed cycle may consist of only one vertex with a directed 
loop. Such a graph would represent a self-duplicating component, 
perhaps a ‘‘living molecule.’’ But nothing more could be said 
about any relational properties of such a graph. From a systemic 
point of view we must have at least two components, the graph 
being then a directed cycle with two vertices or, as we shall call 
it, a 2-cycle. 

In general, as remarked previously (Rashevsky, 1954), any di- 
rected n-cycle corresponds to a system which cannot appear spon- 
taneously. The relational properties of such a graph and of the 
corresponding system are again trivial. 

We must, however, remember that self-reproduction and impossi- 
bility of spontaneous generation are not the only essential char- 
acteristics of a living organism. A system which would reproduce 
itself only under very strictly specified environmental conditions 
hardly would be able to multiply indefinitely because the environ- 
mental conditions never remain absolutely constant. A living sys- 
tem must retain its ability to reproduce even when environmental 
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conditions vary within a certain range. The wider the admissible 
range of variations, the more stable and more viable the organism. 
Its chances of survival increase with the permissible range of 
variation. 

This adaptability of the organism is closely related to the fact 
that an organism has the property of selecting from the environ- 
ment the necessary building stones for its own components even 
when those building stones are present not in a free state but as 
parts of other environmental units. The relational similarity be- 
tween the logical process of selection and the biological process 
of digestion has been emphasized and discussed elsewhere 
(Rashevsky, 1956a). Digestion of food, in the broadest sense, is 
the prerequisite to any further anabolic processes, including that 
of reproduction, Though usually not designated as such, the cy- 
tolysis produced by a virus is essentially a variety of digestive 
process. 

Hence one of the components of the organism must have the 
property of ‘‘extracting’? from the environment the necessary 
building stones for itself and for the other components. Let us 
designate those building stones by Ar Ae ia A. The differ- 
ent components of the organism, in particular the different struc- 
tures involved, which we shall designate by S, are built of different 
combinations of the 4,’s. But the A,’s do not enter spontaneously 
into such combinations which form the different components of the 
organism. Otherwise, since the A,’s may be present occasionally 
in a free state, an organism could be formed Spontaneously at 
such occasions. The presence of an appropriate selective enzyme 
E,, for the formation of the &th component is necessary. These 
enzymes cannot be present in the environment, because otherwise 
again spontaneous generation would be possible. Hence at least 
some of the enzymes E,’s must themselves be manufactured from 
the 4,’s. We may have, for example, the following scheme: 

The environment, Env. plus an enzyme £, which exists only 
within the organism, yield the necessary A.’s. The A,’s plus an- 
other enzyme E, which also is present only within the organism, 
yields the structural element S of which the organism is made, 
This structural element must then yield directly or indirectly both 
FE, and E,. Each process is made through an enzyme which may 
be present in the environment. We shall denote them correspond- 
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ingly by Ey and &’. Thus we have a scheme: 


Env. +E, — A,’s, 


apt E > S, 
(1) 


Shih, 


Se ena Oe 


If we represent the scheme (1) by an oriented graph we see that 
this graph is not a simple directed cycle (Figure 1). The vertices 
which correspond to different parts of the environment are all of 
degree 1. The vertices which correspond to the components which 
receive inputs from the environment are of degree 3. But the com- 
ponent S which represents the ‘‘substance’’ of which the organism 
is built, is represented by a vertex of degree 4. Even if we con- 
sider a simpler scheme, consisting only of enzymes present in the 


FIGURE 1 
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organism and omit the ‘‘structural’’? component S, 
Env. +E =) A's; 
A;'s +E, —> Es (2) 
ASE eee OP 


we find that while the vertex Env is of degree 1, E, and E, are 
each of degree 3, EF, is of degree 4, and 4,—of degree 5. (Fig. 2). 

Without discussing it any further we wish to pose here the 
following problem: What are the general relations between the 
number of vertices, their degrees, the number of organismic com- 
ponents and the number of environmental parts in a graph which 
represents a self-reproducing organism which is incapable of 
Spontaneous generation? 

In other words what should be made is a study of topological 
properties of graphs which represent systems that possess the 
minimum number of attributes which make them classified as liv- 
ing. Such a study may eventually enable us to establish the 
structure of the primordial graph from purely theoretical consider- 
ations. If the solution of the problem is not unique, we shall be 
guided by considerations of relative simplicity. Once the struc- 
ture of the primordial is known, we can derive the graphs of higher 
organisms, for example, by using the postulates discussed pre- 
viously (Rashevsky, 1958). 

The solution of the problem may perhaps be best approached by 
using Robert Rosen’s (1958b) new concept of an ‘‘inverted 
graph,’’ based on the theory of categories. 


FIGURE 2 
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Robert Rosen (personal communication) points out that the sys- 
tems (1) and (2) may be considered as special cases of his (M,R)- 
systems for which he derived a number of biologically interesting 
theorems (Rosen, 1958a). In fact we may consider EF, and E, in 
(1) as two M-components. The A,’s are outputs of E, and inputs of 
E,- We may consider § as the output of E,, and input into EY and 
Es. The first then plays the role of the R-component which corre- 
sponds to E,, the second—the role of the R-component which cor- 
responds to E,, because E, and E, are the outputs of Ej and Ey 
respectively. A similar consideration holds for the system (2). 

It is true that in Rosen’s original presentation (Rosen, 1958a) 
there seems to be a tacit implication that the R-components are 
part of the organism and not of the environments as are E, and Ey. 
However, the above mentioned formal analogy may be significant 
and may form a basis for the study outlined above. 

We may consider two systems the graphs of which are homeo- 
morph, but such that one produces an enzyme eek the other an 
enzyme BS. Let Eby be able to extract the A,’s from a different 
environmental structure than E{*, In biological terms, while one 
enzyme digests only one type of food, the other digests only an- 
other type. The two systems, taken together, form a third, more 
complex system which can use either of the two foods or both. It 
is therefore more adaptable and has a better chance of survival. 
The new system maps onto each of.the old ones two-to-one. The 
two original systems may have some identical components. In that 
case the mapping of the combined system onto each original one 
is not everywhere two-to-one. Combinations of many such systems 
can be formed, and they will satisfy the law of biological mapping 
(Rashevsky, 1959b). The above considerations even suggest that 
the principle of biological mapping may be a consequence of the 
property of self-reproduction and adaptability of biological systems. 

System I may produce substances which are used by system II 
and vice versa. Such systems may also combine and the combined 
system maps on each original one. We may have here the begin- 
ning of specialization and division of biological functions. 

All the above considerations do not, however, answer the ques- 
tion as to how the first organism was formed. A. I. Oparin (1957) 
speaks of the ‘‘gradual evolution’? of the property of self-repro- 
duction but does not give any hints of the specific mechanism by 
means of which a system which is incapable of spontaneous gen- 
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eration could have been initially formed. We shall now discuss 
this point. 

If, for example, in the system (1) the component E, can be 
formed spontaneously, then the whole system can also be formed 
spontaneously. Once formed it could multiply either by spontane- 
ous formation or by self-reproduction according to scheme (1). lf; 
as Oparin suggests, the self-reproduction is a much more rapid 
and biochemically more efficient process, then the self-reproduc- 
tion will gradually predominate and will be much more noticeable 
than the process of spontaneous generation. Still the latter would 
be occasionally present. It is of course possible that the self- 
reproducing system being more adaptable will continue self-repro- 
duction under changed environmental conditions under which spon- 
taneous formation becomes impossible. Still, at least once in a 
while conditions necessary for spontaneous generation should oc- 
cur and the latter should be observed. All this applies not only to 
system (1) but to any other type of such systems. The fact that 
this seems never to happen makes the following conclusion logi- 
cally unavoidable: Once the system (1), or another of a similar 
type is formed, it inhibits the spontaneous formation of E,, or of 
any of its other components. Life, once formed, inhibits sponta- 
neous generation. 

It may be thought that even so spontaneous generation would be 
observable now in such regions of the earth which have not been 
subjected to the effect of organisms. But inasmuch as a sponta- 
neous generation probably occurred from organic material which 
was formed abiogenically (Oparin, Urey) and inasmuch as at 
present practically all organic material on earth has been formed 
biogenically, we should not expect to observe under present con- 
ditions any spontaneous generation. From this point of view the 
recent warnings that in the interest of the study of the origin of 
life we should take special precautions to avoid contamination of 
the moon with terrestrial organic material in projected moon 
probes become of particular importance. 

Though the above mentioned ‘‘antibiogenetic effect of’ life’’ 
seems to follow from an empirical fact, there may be an internal 
logical reason for its existence. Robert Rosen (personal commun- 
nication) remarks that the possibility of representing an organism 
by an oriented graph means essentially that those parts of an or- 
ganism the representations of which do not contain cycles are 
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basically partially ordered sets. In our studies we have thus far 
considered a partial ordering in time. But spatially an organism 
also can be considered as a partially ordered set. The partial 
ordering in time could be upset if some cell-type which represents 
an element of the partially ordered set could be produced sponta- 
neously, so to say ‘‘out of turn.’’ Only if it is produced ina 
proper order from a preceding element is a partial ordering as- 
sured. The impossibility of spontaneous generation may be a con- 
dition for the possibility of development of highly complex or- 
ganisms. 


The author is indebted to Mr. Robert Rosen for a discussion of 
this paper. 

This work was aided by grant RG-5181 of the United States 
Public Health Service. 
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FEEDBACK IN PHYSIOLOGICAL SYSTEMS: AN APPLICATION 
OF FEEDBACK ANALYSIS AND STOCHASTIC MODELS 
TO NEUROPHYSIOLOGY 
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In the human, the antagonistic, extensor-flexor system of the leg is an 
example of a common type of neurophysiological feedback system. After 
a brief introduction to the neuroanatomy and physiology of this feedback 
system, the paper formulates transfer functions from temporal response 
data available in the literature. A feedback stability analysis, based on 
the extension of Nyquist’s stability criteria to multiple-loop systems and 
utilizing flow-graph techniques, demonstrates the stable behavior of the 
system. Expressions are given relating the sensitivity of the system to 
variations in muscle response and Golgi tendon organ (tension receptor) 
response. By considering the events taking place at synapses and end- 
plates during ‘‘isometric tension-small knee angle excursion’’ conditions 
as stationary stochastic processes, an external ‘‘noise’’ input to the 
system is given, whose spectrum is derived from the statistics of a shot- 
process representation of these events. The paper concludes with some 
correlations between the analytical results and clinical syndromes. 


I. Introduction. Feedback is a common occurrence in physiologi- 
cal systems. In man, its macroscopic manifestations range from 
the ‘‘one-shot’’? feedback of reflex response to noxious stimuli, to 
the performance of skilful motor activities requiring the integration 
of a multitude of internal and external stimuli and responses. Mi- 
croscopically, the central nervous system is rife with internuncial 
pathways, centrifugal and lemniscal tracts, recurrent collaterals, 
transcortical and callosal fibers, and so on. The pathways can be 
completely internal or have an external link; with either electrical, 
chemical, or mechanical modalities of operation, or various combi- 
nations of these. Teleologically, the purpose of the various feed- 
back mechanisms ranges from the great homeostatic regulation of 
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the internal environment and the preservation of the bare spark of 
life to the fineness of mental and motor control which enables 
virtuosi to enrich society culturally and aesthetically. 

The various feedback mechanisms associated with motor activity 
are interesting when viewed with the possibility of applying some 
of the available analytical techniques. One of the most externally 
simple feedbacks involving motor activity is the bisynaptic reflex 
arc. This is usually a highly transient feedback situation enabling 
the organism to escape noxious stimuli and their presumed harmful 
effects. The stimulus impinges upon sensory receptors which are 
directly linked via afferent nerve fibers, a synaptic junction and 
an interneuron in the spinal cord, to efferent fibers and an effector, 
such as a muscle. The muscle reacts suitably, removing the or- 
ganism or a part of the organism from the stimulus, and breaks the 
feedback loop. In such a response, fineness of motion is not re- 
quired, and the internal mechanics of the system would be rela- 
tively simple, were it not for the problem of antagonistic muscula- 
ture and reciprocal innervation. In addition, there is the need to 
limit the motion and travel of the member so as to prevent damage 
to joints, muscles, tendons, and so forth. 

Antagonistic musculature is necessary for positive motion and 
control of position and rate of movement of the limb common to the 
agonist-antagonist muscles. Without it, all movements would be 
quite jerky, and we would be most dependent upon gravity effects 
and muscle elasticity for the return of a limb to its original posi- 
tion. The reciprocal innervation is intimately connected with an- 
tagonistic muscles, in that the proper sequences of control of the 
contractions and relaxations of the antagonists is necessary to 
achieve the desired result. Flexion of a limb about a joint requires 
relaxation of the associated extensor muscles as well as contrac- 
tion of the indicated flexor muscles. In more complicated motor ac- 
tivity, such as posture maintenance and stereotaxic motion, it is 
frequently necessary to have simultaneous contractions or relaxa- 
tions in antagonists, such as occurs in the ‘‘locking’’ of a joint or 
maintaining a fixed angle between articulated limbs. The same 
sort of reciprocally innervated effects can occur between entirely 
separate limbs, such as the two legs during walking. Pressure on 
the ball of the foot making contact with the ground causes the lock- 
ing of the knee joint in the corresponding, or ipsilateral, leg and 
unlocking and flexion in the opposite, or contralateral, leg. These 
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latter examples are obviously not simple feedback paths such as 
the bisynaptic reflex arcs described previously and are character- 
ized by complex multiple-loop configurations and a high degree of 
integrated control from the more central portions of the nervous 
system. 


Ul. Neurophysiology of the flexor-extensor system. As an ex- 
ample of a type of feedback system which is neurophysiologically 
common, consider the antagonistic, extensor-flexor system of the 
leg. This system is of major importance in posture maintenance 
and locomotion, and it participates in monosynaptic reflex activity 
such as that found in the stretch, or myotatic, reflex. Figure 1 is 
a simplified neuroanatomical representation of the existing system. 
The simplification resides mainly in the uniqueness of the struc- 
tures. Only one large extensor muscle, the quadriceps, and simi- 
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FIGURE 1. Innervation of the knee flexor-extensor system. A simpli- 
fied neuroanatomical representation of the actual system, in which one 
extensor muscle, the quadriceps, and one flexor muscle, the semitendino- 
sus, are shown. Corresponding muscle spindles, tendon organs, afferent 
and efferent nerves, and representative spinal segments are also indicated. 
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larly only one large flexor muscle, the semitendinosus, is shown. 
Correspondingly, only single muscle spindles or ‘‘stretch recep- 
tors,’’? Golgi tendon organs, and associated afferent and efferent 
nerves are indicated, with a representation of several spinal seg- 
mental levels compressed down to one. 

The alpha and gamma are neurophysiological designations re- 
ferring to the relative diameters of the particular nerve fibers, the 
alpha_fibers are the largest group, ranging from 10-20 p in diameter 
and the gamma fibers somewhat smaller, with a range of 4-10 ux. 
The designations are not too definitive, and the associated fiber 
diameters tend to overlap. Since nerve conduction is an active 
process which depends upon ion transport across membranes, there 
are significant time delays involved in the propagation of nerve im- 
pulses, which are inversely related to the fiber diameter (Lloyd, 
1955). The occurrence of the various fibers in a large- and a small- 
diameter group is functionally important, and the innervation of 
feedback systems of the flexor-extensor type is frequently referred 
to as the alpha-gamma system. The larger-diameter alpha fibers, 
being capable of faster impulse propagation, can transmit higher 
rates of information than can the slower conducting gamma fibers, 
but are relatively bulky. For these reasons their occurrence in the 
nervous system is associated with teleologically important func- 
tions such as those participated in by the system under consideration. 

The muscle actually consists of a great number of contiguous, 
individual fibers which are capable of contraction when electri- 
cally excited. Since not all the muscle fibers in a gross muscle 
are innervated by a single nerve, it is possible to have only a por- 
tion of the fibers in contraction or relaxation with respect to other 
fibers and groups of fibers in the muscle. The individual fibers 
can be only either fully contracted or fully relaxed, their tension in 
contraction a function of temperature and fatigue. A graded re- 
sponse from the gross muscle is thus the result of temporal and 
spatial summation of the ‘‘all-or-none’’ action of the individual 
fibers. The single motor nerve fiber which innervates a group of 
muscle fibers does so by means of the muscle endplate. There is 
no continuum of neuroarchitecture between the nerve and muscle 
fibers, but only a contiguity in the form of the endplate membrane. 
An efferent nerve impulse liberates a minute quantity of chemical 
transmitter substance which diffuses across the membrane and ex- 
cites the postmembranic muscle components into electrical activity. 
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The transmitter substance, acetylcholine, is extremely short-lived, 
being removed from the site almost immediately after its iibersnen 
by the enzyme cholinesterase. 

Associated with each group of muscle fibers in the main muscle are 
the muscle spindles, or stretch receptors. These consist of a group of 
interfusal muscle fibers which are independently innervated, and a 

.neural bag containing tension-sensitive receptors. The spindles 
are effectively in parallel with the main muscle, and length varia- 
tions in the main muscle are reflected as variations in the tension 
applied to the receptors of the neural bag. The receptors generate 
impulses at a rate proportional to the tension applied, and the im- 
pulses are propagated along the small-diameter gamma afferent 
fibers. The tension on the neural bag can also be varied by chang- 
ing the contractile state of the intrafusal muscle fibers. Thus the 
output of the spindle is a function of the relative contraction be- 
tween the main muscle fibers and the intrafusal fibers, and, by 
changing the intrafusal contraction, it is possible to ‘‘bias’’ the 
response, 

In addition to the muscle spindles, there are tension-sensitive 
receptors, known as Golgi tendon organs, located in the muscle 
tendons. These respond equally to changes in tension of the main 
muscle of both a positive and a negative sense, that is, during both 
contraction and extension by the antagonist. The tension threshold 
is somewhat higher for the tendon organs than for the nuclear bag 
endings. 

The motor neuron, located in the ventral horn of the spinal cord 
gray matter, with its associated synaptic endings from other points 
in the nervous system, operates by means of electrical and chemi- 
cal modalities somewhat analogously to the muscle endplates and 
fibers. Impulses arriving at the synapses probably produce the 
characteristic release of a transmitter substance, which causes the 
postmembrane portions of the neuron cell body to become electri- 
cally active. It is not definitely known whether or not the trans- 
mitter is acetylcholine. The integrated effect of the numerous post- 
synaptic potentials thus produced causes the discharge of an im- 
pulse along the motor neuron axon (efferent fiber) if the threshold 
of the electrically excitable portion of the neuron is reached. The 
process may be compared to a shot process, where the randomly 
distributed arrivals of the postsynaptic potentials are due to the 
spatial distribution of the numerous synaptic endings on the cell 
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body and the random interarrival times of impulses at the various 
synapses. The postsynaptic potentials can be excitatory in nature, 
tending to depolarize the electrically active portions of the neuron, 
or inhibitory, with a hyperpolarizing effect on the neuron, depending 
upon the particular synaptic ending involved (Eccles, 1957). 


Ill. Feedback analysis. Over its range of possible behavior, the 
system under discussion is definitely non-linear, and a meaningful 
analysis would be quite difficult. However, some insight may be 
gained by making certain simplifying assumptions. Under the con- 
ditions of nearly isometric tension in the muscles, such as might 
occur during the maintenance of a fixed angle between the thigh 
and foreleg, and where this condition has prevailed for a sufficient 
time to consider the events taking place at synapses and endplates 
as stationary stochastic processes, it seems plausible to apply 
the techniques of linear feedback system analysis. 

With this objective in mind, a block diagram of the feedback sys- 
tem representing the neuroanatomical configuration of Figure 1 has 
been postulated (Fig. 2). Although in a strictly isometric situation 
there would be no length differential between the muscle spindle 
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FIGURE 2. Block diagram of the knee flexor-extensor system as pre- 
sented in Figure 1. Included in each nerve path is a delay element with 
a delay equal to the mean value of a normal distribution of delays ex- 
pected for the particular path. The facilitory or inhibitory nature of the 
inputs to the motor neurons (ventral and gamma horn cells) is indicated 
by a corresponding plus or minus sign. 
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and the main muscle, there actually is some length variation in the 
main muscle due to the elastic properties of the tendons. A spindle 
pathway has therefore been included. In the initial analysis, the 
ventral horn cell (motor neuron) is assumed to act as a linear, at- 
tenuating summing point for the incoming quantities. The neurons 
in the ensemble represented by the single ventral horn cell are 
threshold devices. Under the usual operating conditions, their 
average output is some fraction of their average input. By using 
the ergodic properties of an ensemble with stationary, stochastic 
statistics, the time-average properties of the individual elements 
can be related to ‘‘spatial’’ properties of the ensemble and its 
representative ventral horn cell. In a later section a degree of 
sophistication will be introduced into the representation of the 
ventral horn cell by considering an external ‘‘noise’’ input whose 
spectrum is derived from the statistics of the shot-process repre- 
sentation of synaptic action. 

Figure 3 shows the isometric twitch response, hy (¢), of a muscle 
and a transcendental function which closely approximates it. This 
is the function, 
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FIGURE 3. Muscle and receptor temporal response curves. (a) Typi- 
cal muscle twitch response. (6) A function which closely approximates 
(a). (¢) A function which approximates the response of a tendon organ or 
muscle spindle to a ‘‘suddenly applied stretch.’’ (d) Hypothetical time 
course of the ‘suddenly applied stretch.’’ (a), (6), and (c) are normalized 
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where K,, is the muscle transfer function gain constant. A twitch re- 
sponse is one in which the excitation of the muscle is exceedingly 
short in comparison to the ensuing response, approximating closely the 
impulse response. The corresponding frequency function, M(s), is 


M(s) = Ky E _ Vis aye (=) (2) 


where 


Erfo(X) = ae e-™ du (3) 
x 


is the complementary error function. From (2) we may obtain the 
magnitude as a function of the real frequency, 
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where ,F,(a,6;c) is Kummer’s confluent hypergeometric series. 
The function (4) and its associated phase are given in Figure 4. 
The phase characteristic is one derived from the plot of equation 
(4), treating the muscle as an active source of very large band 
width in cascade with a realizable, minimum-phase shaping function 
not unlike the characterization of vacuum-tube circuits. This is 
based upon the Hilbert transform relationship between the real and 
imaginary components of a complex function and its extension to 
attenuation-phase relationships in networks (Bode, 1945). 

It is possible to obtain data regarding the response R(t) of a 
Golgi tendon organ or a muscle spindle to a “*suddenly applied 
stretch,’’ and this appears to be fairly closely approximated by 
functions having the form (Granit, 1955), 


R(t)=K\Vt ime (5) 
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FIGURE 4. Magnitude and phase angle of the muscle transfer function 
versus real frequency. M(q@) is the Fourier transform of the approxima- 
tion to the muscle twitch response given in Figure 3(b). The phase angle 
®y (@) is related to the Hilbert transform of M(w). The frequency scale 
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The ‘‘suddenly applied stretch’’ could be treated as an ideal step- 
function of tension or, more reasonably, as a function that rises to 
its final value in some finite time. Further, we know that the ten- 
sile properties of the receptor attachments are the result of the 
averaged effects of a large number of individual tensile properties 
such as might be the situation with a group of muscle fibers. 
Therefore, it might be judicious to assume a ‘“‘stretch’’ function 
related to a normal function of time, 
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Both equations (5) and (6) are shown in Figure 3. The frequency 
response corresponding to the muscle spindle (or Golgi organ) im- 
pulse response, SP(s), will then be 
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and 


| SP(o)| =Kgp (8) 


The Golgi tendon organ response, G7 (wa), differs from (8) only in 
the multiplying constant. As before in the case of the muscle re- 
sponse, the phase characteristic of (8) is derived directly from the 
amplitude characteristic. Figure 5 shows (8) and its associated 
phase as a function of frequency. Earlier in the paper, mention 
was made of a difference in threshold between the muscle spindle 
and the Golgi tendon organ. In a general analysis in which appreci-~ 
able amounts of length and tension variations occurred, the effect 
of the threshold difference would have to be included. This might 
be done by means of a describing function derived from the Fourier 
series representation of the nonlinearity. However, in the present 
analysis dealing with isometric tension, it is assumed that the con- 
tribution of length changes upon the muscle spindles to the events 
at the ventral horn cell are small compared to those due to tension 
changes on the Golgi organs. The threshold differences will there- 
fore be counteracted by the relative importance of the contributions; 
the actual tension change in the spindle neural bag being reduced 
by a factor a over those in the Golgi organ. 
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FIGURE 5. Magnitude and phase angle of the muscle spindle or tendon 
organ transfer function versus real frequency. The spindle transfer func- 
tion, SP(@), is assumed to differ from the tendon organ transfer function, 
GT (q@), only by a constant multiplier (gain constant). 
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The delays incurred in the various portions of the feedback path 
are a function of the path length and conduction velocity. These 
paths are representative of an approximately normal distribution of 
similar paths to the various groups of muscle fibers, stretch and 
tension receptors, and motor neurons that are represented as sin- 
gle entities in the diagram. Table 1 gives some figures based on 
an assumed normal distribution. 


TABLE 1. 
Velocity Mean delay 
Path Length (m.) Fiber diam (1) (mps) (ms) 
M S.D. M SHID) M S.D. 

Femoral n. OLD0eeeOL04 815.0 alpha 1.4) 100-0) 1050 5.0 
(spine- 7.0 gamma 1.0 30.0 7.0 17.0 
quadriceps) 

Tibial div. of Same as above 
sciatic n. 

(spine to 
hamstrings) 

Intersegment. 0.02 0.004 10.0 10 5050 958.0) 074 tom.0 
fibers (spi- 
nal Lo_4 to 
L4-S1) Se 

Recurrent milli- — _— —- — — Negligible 
collaterals meters compared 

to other 
delays 

Recurrent (with synapse, i.e. Renshaw cells) 0.5 to 1.0 
collaterals 


In the analysis the mean delays will be used with an impulse 
response, 
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where § is the Dirac delta-function, and complex frequency re- 


sponse of 
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The real frequency amplitude function is unity for all frequencies, 
and the phase angle is 


©, (@) = — Tmean © (11) 


mean 


The system visualization of Figure 2 has been redrawn as a 
flow-graph in Figure 6, in order to facilitate the stability analysis. 
The node of interest, T (tension), has been split and the supra- 
spinal inflows neglected on the assumption of stationary input 
distributions. The question of interest is, ‘‘What sort of behavior 
can be expected of the tension in response to variations of the 
tension?’’ After some reduction, the resultant single-loop flow- 
graph is given in Figure 7. The effect upon over-all stability of 
poles and zeros of internal-loop gain functions cannot be readily 
determined from an examination of the singularities of the resultant 
single-loop function. An analysis based upon the general stability 
criteria of Bode (1945) must be used. In this, the analysis is be- 
gun with one of the internal loops, with all other parts of the path 
reduced to zero. The net number of poles in the right half of the 


FIGURE 6. Flow-graph representation of the system. (a) Initial flow- 
graph with topography exactly as in the block diagram of Figure 2. (0d) 
Unreduced flow-graph, with rearranged topography. The node of interest 
(tension) has been split and taken as the input and output. 
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FIGURE 7. Flow-graph representation of the system. (a) The flow- 
graph of Figure 6(0) after some reduction. (0b) The same after final re- 
duction into a single-loop form. The relations between the reduced ele- 
ments and the system parameters are also shown. 


complex frequency plane is then determined from the encirclements 
of the critical point by the internal-loop gain versus phase plot. 
This procedure is continued by restoring the nulled portions of the 
complete system one-by-one and reanalyzing the growing active 
portion at each restoration. The ultimate stability of the over-all 
system depends upon there being no net counterclockwise encircle- 
ments of the critical point in adding up the results of each succes- 
sive mapping. 

Applying the method outlined above to the system at hand, a con- 
venient starting point is the element G, of internal loop G,. Here 
G, is a loop itself, with the form 


M-D, MD, 
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The internal-loop gain, G,, is plotted in Figure 8. Under normal 
circumstances G, will never encircle the critical point, since even 
the marginal case where G, goes through —1, 0 requires G,(0) to be 
about 300 db. Normal values of G,(0) are about 50-60 db, being 
composed of 40-50 db for K,,, and 6-8 db for K,7. R. Granit (1958) 
has experimentally determined the corresponding value for the cat 
soleus muscle feedback loop of 30-40 db, which, by dimensional 
scaling, indicates that the values given above for G,(0) are of the 
correct order of magnitude. Thus G, can be considered stable for 
all cases. The next internal loop, G,, is composed of G, and G,, 


G 


6s (13) 
i= G, Gs, 
where 
1 
mle ie (14) 
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FIGURE 8. (a) Magnitude and phase angle of Gy, the loop transmission 
transfer function of reduced element G,, versus real frequency. (6) Polar 
plot (Nyquist diagram) of Gy, showing critical point at (-—1, 0). For the 
case where G, just goes through the critical point, the over-all gain con- 
stant of Gy) must be approximately 300 decibels. 
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The behavior of G,-@, about 1, 0 is of interest this time. If M-GT 
is much larger than 1, as is the case, G, is approximately equal to 
1/GT -D,,, and 


(15) 


Thus G,G, is a constant 1/a for all frequencies, and G, has no 
poles or zeros for a greater than 1, which is exactly the situation 
that exists. The marginal case when a equals 1 cannot occur under 
the assumptions of isometric tension with a slight allowance for 
tendon stretch, since this would require the intrafusal fibers to be- 
come the tension member, i.e., take over the function of the main 
muscle. This is an obviously unphysical situation. 

Since G, and G, are in cascade with G,, the effects of all can be 
examined concurrently. The amplitude of any delay member, D, is 
one for all frequencies, with a phase angle that is a linear function 
of frequency. Both D,, (the recurrent collateral) and D,., (the 
intersegmental fiber) are elements of internal loops which include 
the linear, attenuating summing-point representation of the ventral 
horn cell. D,. is a part of G, and Dj, a part of G,. Since the 
amplitude of the delay elements is 1, the stability behavior of G, 
and G, is dependent upon the properties of the summing point. The 
attenuation property insures that G, and G, have no poles in the 
right half of the complex frequency plane, i.e., that the magnitude 
of Dae or Lip in cascade with the attenuation component of the 
summing point is strictly less than 1. 

G, is also in cascade with G,, G,, and G,, and, by examining the 
behavior G,G, about the critical point of -1, 0, the stability of the 
system can be determined. Since 


Ge G, 
B1-G,4, 1-G,6@,’ 
(16) 
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Using the result of (15) and the previous analysis of the behavior 
of (1~-G, G,), 


(17) 


which is always greater than ~—1, 0 for a greater than zero. There- 
fore, under the assumptions given, the system is definitely stable. 
There are, however, clinical syndromes in which the extensor-flexor 
system has an oscillatory behavior. These will be discussed in 
the conclusion. 

The sensitivity of the system to variations of a particular portion 
can be determined by flow-graph methods. The original flow-graph 
is put in the so-called standard form: 


Tate = Tae (18) 
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where T is the over-all transmission, Ty the forward transmission, 
T, the direct transmission, T, loop transmission, and 1 the return 
difference. Sensitivity, 
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if T is much greater than T,, the sensitivity is equal to 1/2. Using 
(19), the sensitivity of the system to variations in muscle response 
is 


1 


Sy 
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where G, is that defined in (13). Similarly, the sensitivity of the 
system to variations in Golgi tendon organ response is 


(21) 


See. GT. DG. eG 
(1 + es | 


IV. Shot-process model of the motor neuron. Suppose that the 
effect of an input at a synapse on the motor neuron is to produce a 
postsynaptic potential component F (¢) in the motor neuron at time 
zero. If we assume that the postsynaptic potential (PSP) com- 
ponents add linearly in the cell, the total effect at time ¢ due to all 
the PSP components is 


+co 
PSP (t) = Ne F(t-t,)> (22) 


k =—co 


where the kth PSP component arrives at ¢, and the series is as- 
sumed to converge. Campbell’s theorem states that the average 


value of PSP (¢) is 


PSP(t) =v i ” F(t)dt, (23) 


oo 


and the mean square value of the fluctuation about this average is 


[PSP (t)- PSP (@)]? = vf F?(t)dt, (24) 


—=Oo 
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where v is the average number of PSP components arriving per 
second. By the method of characteristic functions, which will not 
be discussed in detail, it is possible to compute the probability 
density function related to (22) and (23), 


7 


f(psp)= = f oxp (~2-PSP us y f [et4F ~1)dt) du. (25) 


While this is of some interest, its utility lies in its role in de- 
termining the autocorrelation function of PSP(¢) and the related 
frequency spectrum. 

The average spectrum for random arrival times, ¢,, is 


N(w) =2v|S(o)|? + 2PSP (2)? 8(e), (26) 


where v is the Poisson parameter associated with an assumed 
Poisson distribution of arrival times. A likely form of F(¢) is the 
function 


F(t) =K,te"™, (27) 


with 


S(o) -[ F(t) e7##¢ dt (28) 


defined as the Fourier transform of the PSP component. 
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where N(q@) is the spectrum corresponding to the ‘‘randomness’’ 
occurring in the synaptic processes of motor neuron function. 

One way of envisioning this randomness is as an external ‘‘noise’’ 
source feeding into the idealized system of Figure 2 at a point 
following the ventral horn cell. The effect of the random noise 
source on the output variable is easily seen by constructing the 
system flow-graph with the noise source as an input and the de- 
sired variable as an output. The output spectrum of the noise is 
the product of the noise input spectrum and the squared magnitude 
of the corresponding transfer function. 

For the case where the output variable is tension, the squared 
magnitude of the noise-to-tension transfer function is 


a2(1 + 67272 _ 9 e-¥24 cog V2o] 
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(32) 


and the effect of the random process on the tension is the product 
of (31) and (32), 
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V. Conclusion. In Section III the analysis indicated that the 
antagonistic flexor-extensor system is stable under the conditions 
of physically realizable neurophysiological components—that is, 
under real-life conditions. Another important assumption was made, 
that the supraspinal inputs from the higher levels of the central 
nervous system were quasi-stationary time functions for any steady- 
state tension and, as such, did not affect the stability analysis. 
Under normal circumstances, these assumptions give results which 
correlate well with actuality. ‘‘Normal’’ refers to the condition of 
the supraspinal inputs for isometric tension. 

For instance, the stretch or myotatic reflex is an active attempt 
at maintaining the steady-state tension in a flexor-extensor system 
in response to an external perturbation. Spasticity is a pathologi- 
cal condition which is somewhat like an exaggerated stretch re- 
flex. In this, abnormal but still stationary supraspinal input causes 
hypertonicity in the musculature—the latter a result of the gen- 
erally higher average level of depolarization in the ventral horn 
cell ensemble resulting from the abnormal input. In terms of the 
analysis, there is less attenuation in the representative ventral 
horn cell, and consequently the poles of G, and G, are situated 
closer to the region of instability, the right half of the complex 
frequency plane—sometimes so close that a perturbation which in 
a normal system elicits the usual active resistance of the stretch 
reflex will cause a damped oscillatory behavior (clonus). Hy- 
potonicity, paralysis, and other forms of reduced response occur in 
pathological situations where ventral horn cell attenuation is in- 
creased, such as amyotrophia and poliomyelitis. In myasthenia 
gravis, there is a loss of muscle endplate function which corre- 
sponds to a reduction of the muscle transfer function gain constant, 
kK, in the analysis. An analytical statement of the sensitivity of 
the over-all system to variations in muscle response was given in 
equation (20). 

A number of types of oscillatory movements are not due to any 
instability of a system such as was described here. Parkinsonism, 
athetoid and choreiform movements, and the like are the result of 
lesions in the extrapyramidal system and a ‘‘release’’ from inhibi- 
tion of spontaneous, rhythmical movements initiated by the cortex. 

In summary, this paper has been an attempt to apply the methods 
of feedback analysis and stochastic models to one of the important 
physiological feedback systems. The extension of Nyquist’s sta- 
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bility criteria to multiple-loop systems, flow-graphs, and stochastic 
shot-processes were some of the techniques utilized. 


SYMBOLS 


a—effective reduction in tension on spindle for ‘‘quasi-isometric’’ 
conditions 

A, AT—subscripts denoting agonist and antagonist, respectively 

a—designation of a certain range of nerve fiber diameters 

D (t)—impulse response of a pure delay element 

D(s)—transfer function of delay element as a function of complex 
frequency (i.e., ratio of output to input) 

D,—transfer function of « fiber 

D,, —transfer function of y fiber 

Disp—transfer function of intersegmental fiber 

Drco—transfer function of recurrent collateral fiber 

{)—return difference (1 - T) 

§(t)—Dirac delta-function of time 

db—decibels 

Erf—error function (see eq. (6) ) 

Erfe—complementary error function (see eq. (3)) 

F (t)—component of postsynaptic potential, as a function of time 

,F, (@, 6;c)—Kummer’s confluent hypergeometric series 

f(PSP)—probability density function related to PSP (¢) 

G—transfer function of an element in the multiple-loop system 
visualization 

GT(s)—Golgi tendon organ transfer function, or frequency re- 
sponse, for complex frequencies 

GT(@)—Golgi tendon organ transfer function, or frequency re- 
sponse, for real frequencies 

y—designation of a certain range of nerve fiber diameters 

hy (t)—impulse response of muscle (including endplate) 

kK —arbitrary constant 

Ker —Golgi organ transfer function gain constant 

Ky —muscle transfer function gain constant 

K "—amplitude constant of F (¢) 

Kgp—muscle spindle transfer function gain constant 

L ,-S,—denotes location in spinal cord, i.e., lumbar, sacral, etc. 

Hiettndecle transfer function, as a faioiion of complex fre- 
quency; also muscle frequency response 
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M(w)—muscle transfer function, or frequency response, as a func- 
tion of real frequency 

m—meters 

mps—meters per second 

ms—milliseconds 

p—microns 

N (@)—average spectrum for random arrival times ¢, of F (¢) 

v—average number of PSP components arriving per second 

@—real frequency component of complex frequency s 

PSP (t)—postsynaptic potential, as a function of time 

PSP (t)—average value of PSP (?) 

®, (w)—phase shift in delay element as a function of real frequency 

Rk (t)—temporal response of Golgi organ, or muscle spindle, to 
error-function of time ‘‘stretch”’ 

S—sensitivity of system to changes in an element @ 

S(@)—Fourier transform of F (¢) 

SP (s)—muscle spindle transfer function, or frequency response, as 
a function of complex frequency 

SP(@)—muscle spindle transfer function as a function of real 
frequency 

s—-complex frequency (s = ¢ + ja) 

T—tension 

T—over-all transfer function of a flow-graph in standard form 

T,—forward transmission of a flow-graph in standard form 

T, —direct transmission of a flow-graph in standard form 

T,—loop transmission of a flow-graph in standard form 

t—time 

7—delay time 

6—element under scrutiny, in standard form of flow-graph 
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In this continuation of a previous report it is shown how the Volterra 
population dynamics, which underlies the statistical theory, can be based 
on @ variational principle; how the dynamics can be generalized as re- 
gards both the behavior of total populations and migration phenomena; 
and how many directly observable data, such as amplitudes and fre- 
quencies of oscillation of a population, fit into the statistical theory and 
can test it. Such a test is carried out in some detail using the fox-catch 
data of Elton, with a clear indication that the theory is capable of com- 
prehending the major statistical properties of population-time curves. A 
final section sketches an extension of the theory to cover secular varia- 
tions of external conditions such as temperature of the environment. 


1. Introduction. In a previous report (Kerner, 1957) it was shown 
that the classical Volterra equations, which provide a general, if 
somewhat oversimplified, dynamics of interacting species, lend 
themselves to a statistical, or thermodynamic, type of analysis 
comparable to that used for physical systems having a great many 
degrees of freedom (when the number of species is too large to 
allow any detailed description of the individual population varia- 
tions). The analysis led to a number of experimentally verifiable 
or deniable consequences and contained an unforced statement of 
an important proposition in mathematical ecology, founded on ob- 
servation in the field, due to Corbet, Fisher, and Williams. 

The present note aims to amplify the prior one. First, some 
general remarks on the Volterra equations will be made, and a pro- 
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posal advanced to remedy the serious omission, in the Volterra 
theory and in other demographic theories, of an account of popula- 
tion variations in space as well as time. Then it will be shown 
how the numerical parameters entering the Volterra equations— 
which play the role of ‘‘microscopic’’ parameters in the thermo- 
dynamic view, much as atomic masses are the microscopic con- 
stants of gas theory—may be fairly directly observed within the 
large biological ensemble and not through the practically impos- 
sible observation of isolated pairs of interacting species. Next 
comes a demonstration of how some types of statistical data on 
selected populations, principally amplitudes and frequencies of 
oscillation of populations, fit into the thermodynamic scheme and 
constitute tests for it. 

Following this is an actual, if somewhat rough, test of several 
theoretical conclusions against ecological field data. 

Throughout, our tool will be the Gibbs canonical ensemble to- 
gether with the assumption of ergodicity: time averages = ensemble 
averages. The main viewpoint is that the oscillating populations 
of some few out of a great many interlocking species are among the 
primary ecologic observables that must be fitted into a general 
scheme. The question is, What do such population data mean? 
What elements in the data are not simply of empirical interest but 
are predictable and co-related? What, after all, is there to be 
learned from a population-time curve looking like an overblown 
picture of random noise? The statistical Volterra mechanics gives 
some clear and comprehensive set of answers a priort. If the an- 
Swers are not correct, one has at least something definite to rebut 
and perhaps a point of departure for better theorizing. 

In the concluding section there will be discussed briefly the 
non-stationary ensembles suited to the study of quasi-static alter- 
ations of external variables (such as physical temperature) affect- 
ing biological associations. 

2. Generalities on the Volterra mechanics. Volterra’s (1931, 
1937) theory of interacting populations is essentially constituted 
by the system of differential equations, 


dN. 1 
rae tah spr Ser? (1) 


Le 
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for the population numbers N of the coupled species, with a, = 
-Q,,- Introducing the stationary state dN /dt = 0 as in the earlier 
report and again calling g, the stationary population levels satisfy- 
ing 


¢B,+ >) aN, = 9, (2) 


we have shown that equation (1) may be cast in the form 


; 0G 
ee = 
with 
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G=t = (er = v.); (7, = q,B,) 


As previously, it will be assumed that equations (2) have unique, 
positive solutions and hence that there is an even number of spe- 
cies in association and that not all « are of like sign. The quantity 
G is a constant of the motion holding the cardinal position in the 
present thermodynamic development that the Hamiltonian holds in 
the analogous physical problem. 

We should like to show that equations (3) are comprehended 
under a variational principle similar to Hamilton’s principle in 
mechanics. This should not be confused with Volterra’s (1937) 
variational principle for equations (1), which hinges on the clumsy 
artifice of writing these equations as second-order equations in the 
variables X, = j Nat. 

Consider the function 
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where the matrix of I’;, is antisymmetric and non-singular (there- 
fore, of even order) but otherwise unspecified for the moment. The 
vanishing of the variation of the time-integral of A, with fixed end- 
points, i.e., 


bo 
af A dt = 0, (4) 
t 


or, explicitly, 


Replacing [’;, by — -T,;, and then the dummy index 7 by the equally 
dummy j, sia: says 


: 0G 
ae LF U; = c= e 
j " 
Let us invert this set of equations, solving for 0: 
-1 
7 = (r ye (5) 


Now if T is antisymmetric, so also is [7?: , fort ie ~ (s, rminor 
of I’), and if rows and columns of this minor are interchanged, one 
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gets a determinant, D, identical with the r, s minor, A, except that 
its elements have indices the reverse of those of A; hence, if these 
indices in D are indeed reversed, one gets (~1)D coming to be the 
same as A, since the order of D is odd. 

Thus equations (5), which are equations (3) with y,, disguised 
as ([_ rae are summarized in the ‘‘least-action’’ type of varia- 
tional principle of equation (4). Not only is f Ad¢ stationary about 
(5) for first-order variations in the w’s; it is also stationary to the 
second order of varied w’s. We have, in fact, calculating the varia- 
tion of the ‘‘action’’ to the second order, 


5 J ade= [SOE Ty di +605 - Pons 
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The first integral on the right vanishes, as above, just if (5) holds. 
The second similarly vanishes upon writing the first member of the 
integrand as 
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which cancels the second member identically. 
Following formally for a moment a classic pattern, the ‘‘momen- 
tum’* conjugate to the ‘‘coordinate”’ v, is 


“slur 
i 
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and the ‘‘Hamiltonian’’ is 


Woes oy [,,%%+G=G. 


While in a purely formal sense G may therefore be looked upon as a 
Hamiltonian, it must be recognised that the whole Hamiltonian 
procedure of dynamics is here nugatory, owing to the fact that, to 
begin with, the ‘“‘Lagrangian’’ equations P, = dA/dv, are of the 
first order. 

Turning now to a different matter, it will have been observed 
that no use has been made of any particular form for the primary 
integral, G, of the system (3). This invites some comments on 
generalizations of the Volterra scheme. 

Let us hold to the aim of preserving a basic conservation law, 
which in any subsequent statistical mechanics would be of funda- 
mental importance. Volterra, in his 1931 monograph, had already 
extended his original equations (1) to embrace the wider type of 
conservation, = w,(N;) = const, for nearly arbitrary &, of which 


Lae Dvr -P 4 (Z-toeF) (6) 


(referring just to the original equations) is but a particular ex- 
ample. 

It is clear, however, from the standpoint of equations (3) that an 
appreciably more general type of conservation than Volterra’s is at 
hand. Thus, adhering only to the antisymmetry of the y —en- 
gendered initially by the idea of reciprocal binary interactions be- 
tween species but now not necessarily limited in meaning by this 
concept—it is seen that G(v,, v, «++, %,) = constant, irrespective 
of both the functional dependence of G@ on the » and, of course, of 
the dependence of the v on the population numbers N. Volterra 
conservation, in this view, is simply the special choice G = 
G,(v,), ¥, = 0;(N,). The y,, here, of course, need not be constants 
but may be functions of the dependent variables v. 

We may note in passing a characteristic time-reversal symmetry 
of equations (3), namely, that if ¢ > (~—¢), the equations are un- 
changed if also the signs of all y,, are reversed, meaning in the 
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original usage (y,, = & ./B.B,) a reversal of the roles of predator 
and prey, thus a kind of biologic analog to the principle of dynamic 
reversibility. The existence of the integral @ reflects this sym- 
metry in a manner akin to energy-conservation’s mirroring the time- 
symmetry of Newton’s mechanics. 

It is well to remember at this point the original purpose of hav- 
ing transferred from N-space to the more abstract v-space: it was 
to make v-space a satisfactory phase space in which Liouville’s 
theorem holds and Gibbs ensembles may be studied. Sufficient and 
necessary for the Gibbs statistical mechanics are a mechanics, 
equations (3), admitting a Liouville theorem and a suitable con- 
stant of the motion, plus great ignorance otherwise. Let us see 
about Liouville’s theorem in the context of the generalization of 
Volterra, the other requirements having been met. For a Gibbs 
ensemble of phase points in v-space, the equation of continuity is 


ey ay. 0 
ey 


p being the density, V the velocity of phase points, reckoned as 


point-functions in the space. Expanding the divergence, 


op 


D 
- +V. Vp+pdivV¥ =~ +p div =0. 


We can have Liouville’s theorem, Dp/Dt = 0, just if div V = 0: 


dy,, OG a°G 
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This may be counted as a condition on the y or G@ or both. It is 
seen to be not strongly restrictive; for example, it suffices that, 
with arbitrary G, y,, (v,, %, «++ 5 %,) be independent of v. and 2,. 
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The position altogether is that the Volterra scheme may be sub- 
stantially broadened, and without harm to statistical mechaniza- 
tion. In illustration of the possible utility of the extra breadth, 


note that the canonical distribution retains the form 


w-G 
é 


p=e 


Previously, G was fixed as in equation (6) and so gave us a den- 
sity-in-phase, in which the components G@_(v,) were strictly separa- 
ted, so that the probability densities for the separate species were 
independent and of the form (returning to N-space) 


-1 +x 
P(n,)dn,~n, Gen: n. (7) 


(n, =N./q,3 % = 7/0; 0 = “‘temperature”’ of the association). Now, 
however, we may reject the independence by writing, for instance, 
G=22g(v,) h,(v,) and v,=v,(N,), or G= 2% G,(v;) and v, = 
v;(N,,~+.-,N,); or we may keep the independence (very convenient 
mathematically), G = = G;(v,) and 2; = v;,(N;), such as to give any 
desired form to the probability law of type (7), according to what 
may be dictated by observational experience. This is to say— 
under the circumstances of observation in the field, where it is 
practically out of the question to discern directly whatever micro- 
mechanics is actuating the population fluctuations but rather where 
probability statements like (7) are in the offing—that one may 
possibly glean a clue to the mechanics from just such statements, 
and having the mechanics is to have a basis for suggestions as to 
further observations and otherwise unforeseen connections among 
them. Naturally, it cannot be expected that the whole mechanics 
is inferrable from the statistics alone; and it is certainly arguable, 
to begin with, that a differential mechanics can hold in any strict 
sense, though differential equations have had a value in the past 
and are apt to be useful for some time to come. The point is that 
the sense need not be so strict in the light of the coarseness of 
the observations; a statistical mechanics stemming from a crude 
mechanics may well be sufficient for the data and for the predic- 
tion of data. ‘‘A tous ces developpements théoriques, quelles con- 
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firmation l’expérience va-t-elle? Les vérifications expérimentales 
sont 4 peine commencées’’; this holds nearly as truly now as it 
did in 1931. 

For the present we shall keep to the original Volterra scheme 
because it contains the result (7) having observational basis; be- 
cause it proceeds from reasonable, if simplified, @ priord considera- 
tions; and because it is relatively tractable mathematically. 

3. Hint to a theory of migration: migration as a generalised dif- 
fusion process. Many demographic theories eschew a space-time 
description of populations, confining their attention principally to 
time-variations of population numbers. We should like here to 
sketch a few ideas which may give some account in general terms 
of migrational, as well as temporal, behavior; the ideas are not all 
novel qualitatively, but perhaps their quantitative exposition here 
may be. This stands as something of a digression from the statis- 
tical-mechanical discussion (as yet) but seems close enough to our 
previous general ecological considerations to warrant presentation 
in this place. 

The basic assumption is that, looking at a number of species 
distributed over a large region of space, it is possible to select 
much smaller regions Ar = ArAyAz containing AN, individuals of 
species i such that AN,/Ar may be sensibly reckoned on an ap- 
propriate scale of length as a scalar point function, the population 
density p;(r, ¢) of the ith species (r denotes generic position #,y,2 
in space). Since the organisms comprising each species are gener- 
ally in motion of one sort or another, we also introduce the current den- 
sity, j;(r,), such that the net number of 7’s per unit time crossing 
the small area AS, having unit normal vector f, is j -m AS. In other 
words, the view is that of the bacteriologist looking at a broth by 
eye or through the low-power microscope or of the aviator scanning 
in some perspective the field life below from his high perch. 

It is barely more than a truism now to write an equation of con- 
tinuity: let V be a large enough volume bounded by the surface S, 
and FP; the net rate of increase, per unit volume and per unit time, 
of species ¢ due to its self-propagation and interaction with other 
species. Then (/ denoting the outward normal to S$) 


Ss ee Oe fica dS + [ ear, (8) 
ge s Vv 
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or 
(net rate of increase in = (net rate of + (natural rate of 
V of species 2) immigration increase within 
into V through  V, exclusive of 

S) immigration). 


Upon using Gauss’s law for the surface integral, one has 


Pi 


; + div j,=F,. 


L 


This equation is devoid of real contents until R and j are speci- 
fied. For this, explicitly biological hypotheses must be intro- 
duced. As regards FP, to within the validity of the Volterra mech- 
anics, equations (1), we may, for example, place 


1 
ee ee O55 Pj Pir (9) 
ae, 


after recognizing that, properly speaking, the N’s in (1) are really 
densities (spatially uniform in the case of (1)) rather than total 
population numbers. 

Regarding j, the question is, What may be the nature of the 
motions of different organisms? We call to mind here the remark- 
able quantitative observations of Przibram (1913) on the motions of 
infusorians. Przibram’s verification of the Einstein law—mean 
Square displacement ~ time—and his diagram of a paramecium 
track, so strikingly comparable to that of a Brownian particle, con- 
vince one that the motions of individuals are effectively random 
walks; whence the immediate conclusion that a host of the organ- 
isms will behave as does a host of molecules in solution: the motion on 
the whole will be a diffusion (Chandrasekhar, 1934), with a current den- 
sity proportional to the gradient of the population density, j = — DVp. 
Here D stands for a diffusion constant characteristic of the partic- 
ular type of organism in a specified physical environment and has 
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the meaning, from the theory of random motion, D = nr?/6 (n = mean 
number of displacements/sec, r? = expected squared displacement 
at any step). We are omitting, of course, consideration of any 
induced or systematic motion due, for instance, to mass motions of 
the medium containing the organisms. 

To what extent is the random walk basically prevalent among 
other species, including higher animal species? It appears that in 
a qualitative way it may be so common as to be the rule rather than 
the exception. Witness Elton (1930): 

“It puzzled at least one naturalist working in Siberia to find out 
how the squirrels all seemed to know where the richest food sup- 
ply was to be found in any particular year, and how they managed 
to concentrate so successfully upon the good trees. The reason is 
simply that each individual is constantly going to and fro, reacting 
by migration to the conditions that it meets, and that these move- 
ments automatically result in a readjustment of the density of 
numbers in different places.’’ It is interesting to compare this re- 
mark with one of Schrédinger’s (1946) relative to molecular dif- 
fusion: 

“That this random walk of the permanganate molecules, the 
same for all of them, should yet produce a regular flow... is at 
first sight perplexing .... If you contemplate [in a schematic of 
a solution of varying concentration] thin slices of approximately 
constant concentration, the permanganate molecules which in a 
given moment are contained in a given slice, will, by their random 
walk, it is true, be carried with equal probability to the right or to 
the left. But precisely in consequence of this, a plane separating 
two neighboring slices will be crossed by more molecules coming 
from the left [higher concentration] than in the opposite direction, 
simply because to the left there are more molecules engaged in 
random walk than there are to the right.”’ 

One does not have to look far in even offhand observations of 
field life to receive the clear impression that the random walk or 
swim or flight, in a colloquial, if not a strictly technical, sense, is 
at bottom a dominant theme, though conditioned certainly by many 
factors. To the aerial observer, the squirrel in a homogeneous 
milieu may be reasonably expected to be like the paramecium to 
the inquiring microscopist. 

We propose, then, the preliminary working hypothesis that self- 
diffusion owing to random-type motions of individuals is a basic 
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driving force for the phenomenon of migration; mathematically, that 
a principal connection between population currents and densities 


is 
i;= - DV p;- (10) 


While self-diffusion has been investigated to some extent, for 
example, by Skellam (1951), it covers only the case where the dif- 
fusional motions of different species are mutually non-interfering. 
Apart from the ‘‘point-interactions’’ between species described in 
equation (9), we must expect a motional type of interaction as well, 
that recognizes the possible bias, say, of the motion of predator 
toward prey and of prey away from predator. The currents and not 
only the densities must, in general, be coupled. Here one is re- 
minded of the analogous situation in many-component diffusion 
theory (Kirkaldy, 1958), which gives us the suggestion that a suit- 
able type of generalization of equation (10) is 


(oe DPV; (11) 
J 


where generally the diffusion coefficients D;, are p-dependent and 
may be explicitly time- and space-dependent as well. The contents 
of equation (11) must be understood, however, in quite a different 
sense from the physical one, where the coefficients D;, are con- 
strained by reciprocity laws of a very particular sort. For a prey 
species, 1, and a predator species, 2, for instance, one would 
write 


i, =-D, Vp, - 4, Ve, 
(12) 
ip = — Dy Vp, +d, Voy, 


or in terms of a diffusional-type drift velocity ~ j, 
(net drift velocity of 1) = (self-diffusion velocity of 1) + 
(an escape velocity of 1 ~ local gradient 
of 2) 
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(net drift velocity of 2) = (self-diffusion velocity of 2) + 


(a chase velocity of 2 ~ local gradient 
Oral), 


Thus a kind of diffusional race between the species in which, su- 
perposed on the self-diffusion of each species, are drifts experi- 
enced via the local gradients of the others—so to speak, each 
recognizing which way the population breeze of some of the others 
is blowing. 

In this way is introduced an element counter to the purely homo- 
genizing effect of simple diffusion, a necessary element in a satis- 
factory theory of migration. 

The proposals above have some appearance of multicomponent 
hydrodynamics, quite mechanistic in character. Conspicuously 
absent is any reckoning of fluctuations of populations in a pre- 
scribed region. The strength of diffusion laws comes from the 
random motions of immense numbers of individuals; when the num- 
bers are less than immense, the fluctuations of numbers, incom- 
prehensible in a diffusion approximation supplying only more or 
less sharp averages, are of primary significance. One may per- 
haps regard the present suggestions as the tentative mechanical 
scheme that hides, or awaits, a thoroughly stochastic description 
of individual organism self-production, of organism encounters and 
interactions, and of organism motility. 

We mention briefly some simple examples. For a single species 
under uniform physical conditions (bacteria or Protozoa in liquid 
culture) we have, by equations (8), (9), and (10), 


The same equation governs the neutron density in a nuclear reactor 
(in the well-known diffusion approximation for the neutron motion); 
it will be seen, in fact, that reactor and culture are the same in es- 
sential aspects, with the exception of the boundary conditions: 
Vp-ni=0 at culture boundary, p =0 (approximately) at reactor 
boundary. If, indeed, at the culture boundary there were escape or 
absorption of organisms (and ¢ > 0), we should have a ‘‘critical ra- 
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dius,’? just as for a reactor, which, if exceeded, would give ex- 
ponential growth and in the contrary case exponential decay. 

By well-known methods, solutions to the differential equation 
suited to a given geometry and initial conditions may be con- 
structed. Under spherical symmetry, for instance (florence flask 
seeded isotropically), one gets the superposition of diffusion 
modes, 


é 
ae -~Dx_t , s 
: a 
pltrt).= ec 4 Cy de) Cages as 
1 


where 2, are the roots 4.49, 7.85, ... , of z, = tan z,, and a denotes 
flask radius, the C’s being fixed by the initial distribution p(r, 0). 
The first mode entirely controls the total population, 


a 4 
Nena (= [4a oles dr = 65(Enat) ot 
0 


Sighting along a thin diametral cylinder of length 2a and cross- 
section AS, as in a turbidity measurement, the enclosed population 
is 


AN =2AS / p(r, t) dr, 
0 


which can reveal something of what are the dominant modes and 
what is the value of the diffusion constant D. 


A premonition of serious mathematical trouble is to be seen if 
the self-increase ep is traded for the Verhulst ep — ap”, 


—=DV*p+ep — dp. 
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The last term presents a difficulty that is only slightly amelio- 
rated by the use of a linearization procedure allowing a discussion 
of just near-equilibrium. 

As a second example, consider the extension of the Volterra 
periodic predator-prey relationship for two species, taking in equa- 
tion (12) constant diffusion coefficients, 


Op, 2 
2 
Sais Cow Py ~ 2, V" pg = Py — My Py Py, 
Op. 2 2 
we Se Py + dy V" py = — €y Py + A Py Pg 


The non-linearity is awkward even in the cases of spatial homo- 
geneity (all V? terms zero) or time stationarity (all d/déterms zero). 
Hence we limit ourselves for the present to small deviations from 
the temporal and spatial equilibrium state (p, = €,/%, Py = €,/O,), 
writing py = €:/%,+Q;, po = €,/%, +Q@_, and thence neglecting 
Q,9 

a 


fe] 
“1_D, V7Q, - a, V?Q, +440, =0 (A, =O, €/%5) 


Ce] 
“2D, V?Q, +d, V?Q, - 240, = 0 (Ag = Og €1/0; )s 


Decoupling these, it is found that Q, or Q, satisfies 


p32 aq = 
7, +D,) V? at (D,D, + d,dy) V4 Q - Ayd, + Agd,)Q = 0, 


wherein is seen a certain similarity to the wave equation or, more 
particularly, the equation of vibrations in elastic solids. The 
mixed term V20Q/dt, involving self-diffusion, in general will make 


for damping. 
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It is interesting to observe, in the case of negligible self-dif- 
fusion (D,, D, ~ 0), when the sole driving forces for migration are 
chase and escape of predator and prey, that there are purely propa- 


gating plane-wave solutions, without dissipation, 
Q = e ilk : r—@t). 


The frequency-wave-number spectrum follows the dispersion law, 


/ r r 
w = Vdd, hk? + + hs 
d, dy 


As all wave numbers k are possible, a wave-packet of limited 
spatial extent may be Fourier-synthesized from wave numbers in 
the vicinity of some given one, k,, and the packet will advance 
with the group velocity 


eS 7 ty ( hy + Ao/dy [k2 +2,/d, 
CL] eee ro Ae POET ae fi 


We have a hint of how literally a group of organisms can move 
about as a unit. 

Including, now, the self-diffusion, we find damped propagation, 
o=%4-72f, 


1 2 
B = 5 (D, +D,) k 
(hence strong damping of just the larger wave numbers), and 


a? = kt ld Bue Ba) : 
4 14 — ae + k* (Ad, +Agd,) +A, AQ- 
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If, then, D, = D,, the dispersion is as previously, and if d,d, > 
(-10,)" vas every wave number is still possible; but if dd, < 
(D, -D )?/4, the coefficient of k4 is negative, and only a finite 
ee of wave numbers and frequencies is possible. 

4, Applications of the canonical ensemble. We return now to the 
context of statistical-mechanical analysis of population based on 
the original Volterra scheme. The principal aim is to exhibit the 
power of the analysis in statistical questions relating to the time 
fluctuations of one species singled for study out of a great many. 
In using the canonical ensemble, we are assuming ‘“‘thermody- 
namic’? equilibrium, as previously described, in the large as- 
sociation encompassing the single species, that is, that the as- 
sociation has been let run a long time compared with any time of 
oscillation of any species. In particular, the use of a stationary 
ensemble is to signify that the single population-time curve is of 
the nature of a stationary time series, two long strips of the curve 
having the same statistical properties. As previously mentioned, 
the Gibbs averages are to be counted as time averages. 

Let us notice first the earlier results, that, with G as in equa- 
tion (6), N_ = q,, and that, from the canonical average of (dG/dv,)? 
and of v,(dG/dv,), 


1 (N-—N)? /N N 
Se eel } oo (13) 
rd H 0] 


which gave the fundamental meaning of the association ‘‘tempera- 
ture’? @ as well as, in the last equality, a testable relationship. 
Thus z is to be construed as directly available observationally. 
The temperature itself, together with G, naturally is defined in the 
Gibbs distribution only to within a scale factor, as in ordinary 
physical thermometry; for its numerical determination any suitable 
reference state of equilibrium must first be assigned arbitrarily a 
numerical temperature 6,, which, together with the totally station- 
ary state, 0=0, of minimal entropy, establishes the size of the 
‘‘degree’’ of temperature. This means that the parameter, = NB, 
must also be of arbitrary scale, and so it is, upon recalling that 
the Volterra ‘‘equivalent numbers’’ pas have no individual signif- 
icance but rather that only the composite opened ee has direct 
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meaning (number of s’s lost or gained per number of r’s gained or 


lost, per binary (7, s) encounter). 
The quantity 


0G. N 
eae eee -1)=4,(71) 
Ov q, 


r 


is plainly of primary importance, so we record its moments (drop- 


ping the r): 


or, since 0?G/dv2 = dG/dv +1, 


D" =(n-1)0D™-1 4 (n-1)r9 D™-?, 
Dia Oe pMere (14) 


D3 =276?, D4 = 676° + 37267, 2... 


The moments of v = log (N/q) may also be computed from the mo- 
ment-generating function, 


ete P(A + D(A + @) a> 


ET fain, 


’ 
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= (s s (2) =1 
ge arc = P(x) -— log a, 
or ae 


a (15) 
v* = p(x) + [p(x) - log a]?,..., 


where 9 denotes the digamma, 9’ the trigamma function. 
As a second step we consider averages involving derivatives of 
the v’s. We have 


since each D by (14) vanishes; also 


ub 
2 == Se _ 
b, = eé6 > Year's fe Ss v, é 6 dv, eco dv, = 0, 


since each ¥, is independent of v, and each v = 0. The combina- 
tion 0,0, must be expected to tell something of the interaction of 
species p and 7, 


¥ Ee ae 
i hy Bee = fo, (a Ysr a é y dv, spo dv, 
0G 
=p, Pra Vp, Oe 
Pp 


Only the pth term in the ». survives, due to D=0. In similar 


s 
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fashion, 


Rfee v 
pppoe t(é,° = 1) Gir yep 
Pp 


These give us the means of observing the ‘‘microscopic’’ para- 
meters y,, of the Volterra theory, 


ee ee ep 
eeT Bt Ip ON. Ip 


or, in terms of the original configuration of parameters entering the 
Volterra equations (1), 


a. N N N N 
Pr. log 2 fr (t= log —., 
B, N, Tp q, oF 


where the numerator may be replaced by the alternatives in (16) and 
the denominator by the alternative for 6/r in (13). 

A basic meaning for the Volterra antisymmetry on the level of 
observable time averages now is 


vo +O 0 0 eo 0. 


dt 


Several other out of many interesting averages may be mentioned 
without proof: 


so) ig. he 0 dam Yir Tie 
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0? = 0 yy Vir 7; ’ 


Pleveal) =~ 0 Der r 
b 


2 


2 
0 2 
: Yir Ti) 
: 


[oun \ = ~ 
" \ao, ) “Yer Op» 
is 


B(eP-1)aQy 767, 


os jones 1)? se 5 


We turn now to the basic problems of gauging something of the 
horizontal spreads of population-time curves. 

The fraction T_/T of a long time interval T spent by a popula- 
tion below the average population level N= q, is the time average 
of 


h(v,) =1, u.< 0 
= 0, v.> 0. 


Using the canonical ensemble, this is (dropping the subscript) 


of he x ) ae 0 Paes i ) oS. 
=- 0 | h(v)e odo= f e av | f e 8 dv 


—0oO 


x co 
= f See eer Sig ih CLE adhe 
) 0 
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after introducing s = ze” as variable of integration. In terms of 
Pearson’s (1951) incomplete gamma-function, 


1( ye 1 eal poate eee 
Uy P pacer G—* -aP aa, 


the mean below-average time is 
Re Za 
erg ie ras 


which varies between 0.5 and 1. It may be emphasized again that 
z is to be considered accessible to observation via equation (12). 
In Figure 1 is plotted T_/T and its complement, the mean above- 
average time 7,/T=1-T./T. At very low association temp- 
eratures, 0 < <7, ise., when the association is not far from the sta- 


FIGURE 1. Mean below-average (7_/7) and above-average (7,/T7) 
times as functions of 2 : 
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tionary state N = q, (all r), the populations spend as much time 
above as below their average levels, oscillating very roughly 
sinusoidally about these levels. For very high temperatures, far 
from the stationary state, the populations spend most of their time 
at below-average levels, oscillating below in long shallow troughs 
and above in short, high peaks. 

To make more precise these views of the amplitudes of oscil- 
lation, we may append to the evaluation of 7,/T, T_/T the cal- 
culation of the separate mean amplitudes A,, A_ of oscillation 
below and above average, averaging separately over the time seg- 
ments when AN > q and when VN < gq: 


Ave fcer—vae [fam f (F-) af fa 


ey rs 
=) (1 - A) (e 1 ae | 


1. 


Cc) 0 


and similarly for A_ if (7,,/T) is replaced by (T_/T). The + on the 
integrals signifies times during which N>gq. Figure 2 gives the 
amplitudes as functions of 2, bearing out the approach to equality 
(and also over all decrease of amplitudes) for descending 0 and the 
sharp increase of A, over A_ for ascending 0. 

A key datum of great practical interest to be comprehended 
within a demographic theory is the frequency of oscillation of a 
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7 aT 2 3. >=4nu Scene t c>-.8 
x 


FIGURE 2. Mean amplitudes of oscillation above (A,) and below (A_) 
average. 


population. It is surprisingly difficult to bring forth this informa- 
tion explicitly, even in the simple case of the periodic two-species 
Volterra interaction, but it turns out to be fairly amenable to treat- 
ment by the canonical ensemble. 

Consider an oscillatory function F(¢). The integral 


- 
/ 5 (F(t) dt, 


0 


where 6 designates the delta-function, will give a contribution 


fe (F(a) (68) dt = 
| F(t.) | 


near a zero ¢=¢, of F, but not otherwise. Therefore, the mean 
frequency of zeroes, w, of F(t) is* 


1) Re 
om f |F’(t)| 8(F(t)) dt. 


0 


*I am indebted to Professor M. Kac for pointing out to me this result 
and its utility in the statistical mechanics. 
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For the mean frequency of zeroes of a population about the level 
N. = 7, Or equivalently about v, = log N/q, = 9, we then have, using 
a Gibbs average for the time average, 


¥ eos 
o,=e® {14 5(v) e 8 dV, +... dv,. 


Since 0. is independent of v,, this breaks up into a simple integral 
over v, and a more difficult one over the rest of phase space dv’, 


Ye G,(0) GG. 
@=e%e 9 pee 6 dv’ 
hae : ee: 
-e6e *|+| =-4+—_|3 |; (17) 


so the computation reduces to that of | % |» 
By a slight extension of the foregoing we also find the mean 
frequency of zeroes of N/g about an arbitrary line N/q = v to be 


N Naot. ©: i 
»,(e ) =e 0 ce od ak d(e7™-v) dv,.+. dv, 
ie 


N N 
= frequency about —=v relative to frequency pees = 1. 
q 


Since a, (v)< 1 if v>1 or vy <1, we find that N/q crosses the 
axis N/q=1 more frequently than any other axis N/g=v. How 


242 EDWARD H. KERNER 


much less is the frequency of oscillation about axes other than 
N = q is to be seen in Figure 3, where contours of a,, = const in 
the z,v plane are shown. It is evident that the frequency of zeroes 
falls off with extreme rapidity about increasingly elevated axes for 
the lower association temperatures (higher x) and that only if 
6>>r7 is the fall-off slow. This is to say, looking also at v< 1, 
that at low temperatures small excursions away from N = g occur 
appreciably more frequently than larger ones (with more small up- 
ward excursions than downward), and contrariwise for high tempera- 
tures. 

The variable w = log N/g is of as much interest as N, and for its 
zeroes-frequency about v = v (v now any positive or negative value) 
is found 


e-x(e — V) [o|. 
r 


Ce 


T(z) 


The rigorous computation of | 2, | unfortunately seems to be out 
of the question, though satisfactory approximations are to be had. 


FIGURE 3. Contours of @., = const in the z, v plane, 
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The magnitude of the difficulty may be seen by writing 


ees! ° 1-cosvé 
|v,| = - ip peace oo dé 


and then 


: co od , : 
ev7e f Bf 87C1 — cos v €) dv’, 


7 


|, 


primes denoting the absence of the variable v. Performing the 


integration over phase space, after placing 3, = = y, 7,(e ‘- 1), 
gives 
| | ey te 
v, = S & ge [ me e (€)] ’ 
where 
sie: 7; ir 
eS Pouce —x. 
LG -~i€Oy,)" i¢r a 


is the characteristic function of v, and by equation (2), = Vie 
ce,» Plainly, v has a complicated distribution, and to find [ 2, | is 
to seek the Fourier transform of a difficultly transformable func- 
tion. From © are visible the semi-invariants of 0, 


k,=0, «,=1! Ay, Ks = 2! ay, 


where y?= Lerey”s y =Xry,?, .+., etc., so that the moments 
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of v, come to be 


9 aa 
v= 0, vu, = Oy", 
2 


0.8 267 y*, v,* =367y? +66°y%, 


a5 220 6 yy? +24 04 y5, ..., 


from which may be determined approximate distribution functions of 
Gram-Charlier or Pearson type. The appearance of the moments 


ys vo ..., raises the interesting side question whether there may 
not be another order of statistics connected with the considerable 
ignorance surrounding the numerous microscopic parameters; 
whether, so to speak, equation 2, (e) = y(r), does not define in some 
way a kind of algebraic statistics on all vectors (e) and trans- 
formations y, admitting only positive vectors (r), over which ca- 
nonical averages would be further averaged. 
For sufficiently small 6 we have 


Q = exp - = [iér, y,, + x; log (1 -7& Ay;,)I 


. . 1 . 
~exp-% ifr y, +2, [-t€Oy,, - 5 GE Oy,) | 


1 ae 
~ exp - = 7 0y, 


and for la] 
ceo cask rere 
| = (*) ov} (2) (52), (18) 


We are using here effectively the fact that at sufficiently small 6 


the distribution of n, = Ye 


%y 
x; —%.R-; ee 


P(e) "i 


u 


P(n) dn, = 
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is sensibly normal, so that 


1 = > ¥;, 5 (0, - 1) = 2 y,, 101+ 


is also essentially normally distributed. 

One has the suggestion here to appeal to the central limit theo- 
rem and to say that, in general, 7, being a superposition of a great 
many independent variates, tends to be normally distributed with- 
out restriction on 6. This may in fact be made quite plausible (if 
not conclusive) in spite of the awkward characteristics of P(n) for 
larger 6. Adopting this suggestion, equation (18) holds generally, 
giving, for the zeroes-frequency of a, 


1 1 
2\zr Be A Bae cng 
(=) Ga oe) 


Hence , is slowly and monotonically varying with 2 (Fig. 4). 
While at @=0 (where the association is completely quiet, ail 
Nt) = q,) the oscillation frequency is zero, as soon as the as- 
sociation is the least bit excited, the frequency jumps to its 
largest value, decreasing thereafter uniformly as 6 is raised. The 
discontinuity in w at zero temperature is, of course, to be expected 
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FIGURE 4. Mean frequency of oscillation @ about N = gq. Ordinate w’ is 
@ of equation (19) with the constant factors omitted. 
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insofar as the perturbation of an initially quiet oscillatory system 
must result at once in a finite frequency, if only an infinitesimal 


amplitude. 
It may be that, hiding behind the smoothness of (19) and covered 


over by the use of the central limit theorem, there is some more or 
less fine structure, a slight waviness, in the variation of » with 0. 

It is interesting to compare the frequency of oscillation found 
above with that which would occur if v(t) were random (Gaussian) 
noise. To this end we first generalize (19) to give the frequency 
of oscillation of » about its average, v (eq. (15)), 


o,(v, - 7%) = 5(v,-%) |%| =5(%,-%) | 2,| 


Seep eens a+ 
z (} (ry2)? Sa ; exp — a[e? 1%* — (9 - log z)]. 
= ! 


Now, if U(¢)=v-—% were a random noise, the mean frequency of 
zeroes would be (Rice, 1944) 


rie ~ (05/09). 


Using canonical averages, 


U0 aleneye eee ets 


U? =(v —0)? =v? — 07 = o’(2); 


whereupon 


es 


| a4 
f =(ry?y ~ (@ 9’ (2)] = 3 


In Figure 5 are plotted f and w(v —%), omitting the common fac- 


wih 
tor (ry")?_ The two results are altogether comparable; the intima- 
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FIGURE 5. Mean frequency of oscillation of v about v =v (curve 
C-C) compared to similar result if » were random noise (curve N-N); 


factor (ry?) has been dropped in both curves. 


tion is that U(¢) is possibly not far different from random noise in 
a general way. 

In passing, it may be mentioned that one may make over the 
population fluctuations into something very apt to be like random 
noise by introducing, instead of v, the function Q(v), defined by 


x 


* exe =v) ree 1 oar dQ, (20) 
T(z) 


aN2n 


that is, by making over the v-distribution into a normal one; this is 
only necessary and not sufficient for Q(t) to be a random noise. 
For large 2(6 < < r) the v-distribution is already normal, 


1 
x Xo X% —_ 1 
e 1 othe 
: ere =e) a . ~zx0" i) zx 


(2) ~ P(x) ew 


(using the Stirling approximation in the last step), so we are guided 


toa = a7F Integrating (20) on the left from 0 to v, and on the right 
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from Q, to Q, Q(v) is then completely defined by 


(Je e", @-1)~ (Vz, 2-1) = £(o V3) : r(o, vz}, 


E denoting the error function, 


1 z 2 
Riavie l= -" dt, 
(2) aa e 


7 0 


The constant Q, is chosen so that Q ranges from (—«, 0) when w 
does, or 


E (0, v3) = I(Vz, 2-1) - = 


5. Comparison with observation. If we settle for a treatment 
grosso modo for ecological field data, some test of a number of 
theoretical conclusions is possible. For the meaningful evalua- 
tion of time-averages, population data extending over a sufficiently 
long time are needed. One of the longest unbroken records, Figure 
6, seems to be the catch of foxes by the Moravian missions in 
Labrador from 1834 to 1925. Elton (1942) has compiled (chap. 13, 
Table 17) and discussed the data. We shall assume simply that 
the fox-catch accurately samples the population, an assumption 
that is not unreasonable, judging from Elton’s elaborate discus- 
sion (chaps. 14 and 15). The catches represent intra-annual 
averages and so, desirably for our purposes, they average over 
seasonal and other smaller-scale fluctuations outside the present 
theoretical framework. Rather than attempt to smooth the data, 
which would be at best somewhat ambiguous and could unduly bias 
the computation of time-averages, we leave them entirely un- 
processed, taking the polygonal line in the upper part of Figure 6 
as the population-time curve. 
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FIGURE 6. Labrador fox-catches for the first 40 years of a 91-year 
period, after Elton (1942). Upper curve gives the catch directly (left- 
hand scale); right-hand scale gives the reduced variable n= N/q. Lower 
curve Shows v(¢). 


In Table 1 we record the directly computed time-averages of 
several quantities, the theoretically consequent values of z, and 
the excursions in the latter due to a +10 per cent alteration of the 
former. 

Aside from the assumption that catch ~ population, errors in the 
reduction of the data are: statistical errors due to the rather small 
time interval of $1 years covering the data; and the errors of the 
assumed polygonal shape of the population curve. Over all, an 
error of 10 per cent in any time average would seem to be a modest 
estimate. The uniformity of x-values, which is the test of the 
theory, is in this light even a little surprising. 


TABLE 1 
Time Average x Excursion 
car 0.595 1.90 0.70-15.0 
ae 445 2.10 1.00-2.75 
A, .729 1d 1.42-2.05 
(n — 1)? 536 1.87 1.69-2.07 
(n—1)logn .489 2.04 1.86-2.27 


log n — 0.258 2.08 2.31-1.91 
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In Figure 7 is shown  |,(v) observationally and theoretically. 
Here the statistical errors are very appreciable, on the order of 
perhaps 50 per cent, since each datum involves the ratio of the 
count of zeroes on two axes, each count having itself a large error 
(maximum count = 43, on axis n=1). The shape of q,,(v) must 
be considered to be satisfactorily accounted for. 

The testing of absolute and not relative zeroes-frequencies 
would be desirable but, as this involves a reckoning of a time 
derivative (v), is scarcely feasible under the ‘‘polygonal’’ as- 
sumption. However, even with this assumption, it is found that 


|| = 1.152 yt 1, w(n—1) = 48/91 yr A giving, according to equa- 
tion (17), «= 1.25, and an excursion induced by a 20 per cent 
error in w/ | 0| of 1.70- 0.84. The polygonal assumption, of 


course, greatly falsifies | | » giving an overestimate; any reason- 


able smoothing decreases || noticeably. For example, assuming 
that » varies linearly over one-year intervals ¢, ¢ +1 between the 


ine) Feaerl (aya 


these v’s are taken right from the data, one gets || = 0.841 yr, 


giving x = 2.16. 
The comparison of theory and observation outlined here is, on 
the whole, more a test of the probability law 7 than of the under- 


x= 15 
X=—"2H 
X= 2.5 


e observed 


FIGURE 7. Comparison of theoretical and observed @,.;(V). Abscissa 


is v in units of 0.16, which correspond to intervals of 50 fox-catches in 
Figure 1. 
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lying Volterra mechanics; a significant exception is the testing of 
®,.; and w, which relies on the basic property of the mechanics, 
that v, is independent of 2,. 

One test does not, in general, make (or unmake) a theory; but we 
have at least the intimation of validity of the Volterra statistical 
mechanics. It is difficult to escape an impression in the present 
case that a Volterra oscillation is effectively, if not literally, at 
hand and that the statistical theory can be a useful tool in the 
interpretation and correlation of ecological field data. Further 
tests would be highly desirable: particularly, for instance, the test- 
ing of d/dt(v,v;,)=0 for two strongly interacting species, that 
hinges so evidently on the basic premise of antisymmetry of y;,- 

6. Non-stationary ensembles. There is one case in which ther- 
modynamically non-equilibrium configurations of macroscopic sys- 
tems can at present be profitably treated within the Gibbs metho- 
dology, that of a sufficiently slow alteration of external variables 
or of the system temperature. We shall adapt the discussion of 
Cox (1950) to a sketch of this case for the macroscopic Volterra 
systems. 

Representative types of external variables of importance in 
biological associations are the physical temperature or radiation 
intensity or oxygen tension in the environment, and the question to 
be answered is, What Gibbs ensemble is appropriate for the de- 
scription of the association when such variables vary and induce 
changes in the microscopic parameters ¢«, B, & (eq. (1))? Addi- 
tionally, the association temperature 6 may be subject to change 
through contacts with a second association at different tempera- 
ture. 

A slight revision of the ‘‘canonical’’ demographic equations (3) 
is first called for so as to make the meaning of phase space inde- 
pendent of changes of the external variahles, X. Let «, B, % be 
(X), B(X), «(X), with X = X(t), and let X, be some fixed reference 
values of the X. Calling 7, = ¢,(X) and q,9 = 7,(Xo) the solutions 
to equations (2) for the indicated values of X, and redefining the 
v’s as v. = log (N,/9,9), the Volterra equations hecome 


a aa 
ee) Peed Beda. ° ~B.9.)= Di Yer Gy? 
x s 
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where y = y(X) and G, no longer a constant of the motion, is 


G=% (B, G0 er oa B, qT, vy, )e 


Imagine, now, that one or more of the X,, or 0, changes so slowly 
that a glance at the association shows it to be nearly in thermo- 
dynamic equilibrium with the momentarily prevailing G(X), giving, 
in the canonical distribution, a close approximation to the true 
density-in-phase. ‘‘Momentarily’’ means a sufficiently long moment 
for a good time-average of some population variable to be taken, 
but so short on the time-scale of the changes in X or 6 that G(X) 
varies only slightly in most parts of phase space. In short, the 
changes envisioned are quasi-static, the system being gently 
propelled through a succession of near-equilibrium states; it is 
only in this case that it makes sense still to speak of a system- 
temperature at all. In seemingly sharp contrast to the physical 
case, the ‘‘moment’’ biologically must be on the order of months or 
years, and the secular changes measured as increments per decade 
or century. 

Following Cox, we suppose that the corrections to the instan- 
taneous canonical distribution are proportional to the adequately 
small velocities of change of the drifting X and @, 


p=e (1. +oAGaigliek camel Mee od 
mip UnetGsael Kei Be Xe Bron (21) 


The A,B are taken to be functions of 6, X, v, but not of A, te 
They are fixed by Liouville’s equation, 


dp ./ Op 
aaa av,’ 
which now reads, after dropping a factor p, throughout, 


1 fap oG\. Y-G-: i Ae 
Daa areata 7 )ardos Bisse 
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On the left, products and derivatives of 6, x may be dropped; on 
the right, the first term vanishes, owing to = v, dG/dv.= 0. In 


what remains we may equate coefficients of the small but arbitrary 
6, X., giving 


—— — —=} = — y —+, 292 
aX, al es av. a 


These are partial differential equations for A,B. Since in equation 
(21) both p and p, represent normalized probability densities, A 
and B , have the restriction A=0, B, = 0, where the bar designates 
an average over the unperturbed senate Poe 

We note a characteristic lack of reciprocity of the Onsager type. 
Define, for convenience, 


5B dB 0B ORs 


= ys 
aE ae dv, A 


and take the p-average of (22), denoted by the double bar, 


F, = ~ —— = average ‘‘force’’ conjugate to X, 
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= eas 5B 5B, 8B, 
ee = =. = a4 xX, ie ae eS ° 
2 =~ 9x o( 222. § sf St 


Also 


oF, Pu es os 5A 
eres Seats pe ee 
06 ax, 
a agey sO Bap Al seveg PBs 
ci] Lee aoe <+ 2. sueeeeee 
aX, aX, 


whereas for physical systems the inequality is an equality. 

The differential equations (22) have been solved in some re- 
stricted cases, but the further elaboration of non-equilibrium awaits 
their more general analysis. 


I am indebted to the National Institutes of Health for a grant 
which supported a segment of this work. 


LITERATURE 


Chandrasekhar, 8. 1943. ‘Stochastic Problems in Physics and Astron- 
omy.’’? Rev. Mod. Phys. 15, 1-89. 

Cox, R. T. 1950. ‘*The Statistical Method of Gibbs in Irreversible 
Change: ” Rev. Mod, Phys. 15, 238-48. 


BIOLOGICAL ASSOCIATIONS 255 


Elton, C. 1930. Animal Ecology and Evolution. Oxford: Clarendon 

Press. 
- 1942. Voles, Mice, and Lemmings. Oxford: Clarendon Press. 

Kerner, E. H. 1957. ‘‘A Statistical Mechanics of Interacting Biological 
Species,’? Bull. Math. Biophysics. 19, 121-45. 

Kirkaldy, J. S. 1958. ‘‘Diffusion in Multicomponent Metallic Systems.’ 
Can, J. Phys. 36, 899-906. 

Pearson, K. 1951. Tables of the Incomplete Gamma-Function (reissue). 
Cambridge: University Press. 

Przibram, K. 1913. ‘*Uber die ungeordnete Bewegung niederer Tiere.’’ 
Pfliigers Archiv, 153, 401-5. 

Rice, S. O. 1944. ‘*Mathematical Analysis of Random Noise.’’ Bell 
System Tech, J. 23, 1-162. 

Schrodinger, E. 1946. What /s Life? Cambridge: University Press. 

Skellam, J. G. 1951. ‘*Random Dispersal in Theoretical Populations.’’ 
Biometrika. 38, 196-218. 

Volterra, V. 1931. Lecons sur la théorie mathématique de la lutte pour 
la vie. Paris: Gauthier-Villars. 

, 1937. ‘*Principes de biologie mathématique.’? Acta Bio- 

theoretica, 3, 1-36. 


RECEIVED 9-15-58 


